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Abstract. The TB-LMTO-ASA method is reviewed and generahzed to an ac- 
curate and robust TB-NMTO minimal-basis method, which solves Schrodinger's 
equation to A'^th order in the energy expansion for an overlapping MT-potential, 
and which may include any degree of downfolding. For N — 1, the simple TB- 
LMTO-ASA formalism is preserved. For a discrete energy mesh, the NMTO basis 
set may be given as: x''^' {^) ~ ^ (^"' terms of kinked partial waves, 
(j){e,r), evaluated on the mesh, eo,---,EN- This basis solves Schrodinger's equa- 
tion for the MT-potential to within an error oc (e — eo) ••■ (e — £n) ■ The Lagrange 
matrix-coefficients, as well as the Hamiltonian and overlap matrices for the 

NMTO set, have simple expressions in terms of energy derivatives on the mesh of 
the Green matrix, defined as the inverse of the screened KKR matrix. The varia- 
tionally determined single-electron energies have errors oc (e — eo)^ ■•• (e — eat)^ . A 
method for obtaining orthonormal NMTO sets is given and several applications are 
presented. 

1 Overview 

Muffin-tin orbitals (MTOs) have been used for a long time in ab initio cal- 
culations of the electronic structure of condensed matter. Over the years, 
several MTO-based methods have been devised and further developed. The 
ultimate aim is to find a generally applicable electronic-structure method 
which is accurate and robust, as well as intelligible. 

In order to be intelligible, such a method must employ a small, single- 
electron basis of atom-centered, short-ranged orbitals. Moreover, the single- 
electron Hamiltonian must have a simple, analytical form, which relates to a 
two-center, orthogonal, tight-binding (TB) Hamiltonian. 

In this sense, the conventional linear muffin-tin-orbitals method in the 
atomic- spheres approximation (LMTO-ASA) is intelligible, because the 
orbital may be expressed as: 

XRL{rR) = (t)RL{YR)+^(j)R'L'{vw){HR>L',RL-£u5wR5L^L)- (1) 

R'L' 

Here, (t>RL (jr) is the solution, ipRi (e^, r^) Yim (ffl) , at a chosen energy, of 
Schrodinger's differential equation inside the atomic sphere at site R for the 
single-particle potential, J2r '^R ('"fl) j assumed to be spherically symmetric 
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inside that sphere. Moreover, r/j = r — R and L = Im. The function ipm (e, r) 
thus satisfies the one-dimensional, radial Schrodinger equation 
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ripRi{e,r). (2) 



In (|f]), (r) are the energy-derivative functions, dtpRi (e, r) /Sel^ Yim (f ) . 
The radial functions, tp and (/?, and also the potential, v, are truncated outside 
their own atomic sphere of radius s, and the matrix, H, is constructed in such 
a way that the LMTO is continuous and differentiable in all space. Equation 
(0) therefore expresses the LMTO at site R and (pseudo) angular momentum 
L as the solution of Schrodinger's equation at that site, with that angular 
momentum, and at the chosen energy, plus a 'smoothing cloud' of energy- 
derivative functions, centered mainly at the neighboring sites, and having 
around these, all possible angular momenta. 

That a set of energy-independent orbitals must have the form (|l|) in order 
to constitute a basis for the solutions if', (r) -with energies in the neighbor- 
hood of £1,- of Schrodinger's equation for the entire system, is intuitively ob- 
vious, because the corresponding linear combinations, J2rl Xrl (t^r) crl,!, 
will be those which locally, inside each atomic sphere and for each angular 
momentum, have the right amount of -provided mainly by the tails of 
the neighboring orbitals- added onto the central orbital's Lp. Since by con- 
struction each ipRi (e, r) is the correct solution, this right amount is of course 
£i ~ Ev In math: since definitions can be made such that the expansion ma- 
trix HniL'^nL is Hermitian, its eigenvectors are the coefficients of the proper 
linear combinations, and its eigenvalues are the energies: 



RL RL 



CRL,: 



^ XRL {yr) CRL,i = ^ ^(pRL (yr) + (£j - e^) (j)RL (Jr) 
RL 

^ 4>RL {£i, yr) CRL,i = I'i (r) . (3) 



RL 



Hence, H is a Ist-order Hamiltonian, delivering energies and wave functions 
with errors proportional to (e^ — £,y)^ , to leading order. 

First-order energies seldom suffice, and in the conventional LMTO-ASA 
method use is made of the variational principle for the Hamiltonian, 

7i ^ -V^ + Y.r'^R^'-i^^ ' (4) 

so that errors of order (e^ — Sy)^ in the basis set merely give rise to errors 
of order (e^ — e^) in the energies. With that approach, the energies and 
eigenvectors are obtained as solutions of the generalized eigenvalue problem: 



X! [(Xi?/L' |W - £1.1 Xrl) ~ (ei - £v) {xr'l' I Xrl)] crl,i = 0, (5) 

RL 
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for all R' L' . If we now insert in (^), we see that the Hamiltonian and 
overlap matrices are expressed in terms of the Ist-order Hamiltonian, H, 
plus two diagonal matrices with the respective elements 

4>RL I (t^RL) = JofRi (r) 'pRi (r) r'^dr, /^rl \ ^rl) = ^Ri {rf r^dr. 



These matrices are diagonal by virtue of the ASA, which approximates inte- 
grals over space by the sum of integrals over atomic spheres. If each partial 
wave is normalized to unity in its sphere: ipm (r)^ r^dr = 1, then (0 | 0) 
is the unit matrix in the ASA, and the Hamiltonian and overlap matrices 
entering (0) take the simple forms: 



(6) 



(Xlx> = 1 + {H -e,){^\^) 1 



2" 



Here and in the following we use a vector-matrix notation according to which, 
for example xrl {yr) and xrl (r^)* are considered components of respec- 
tively a row-vector, x (r) , and a column- vector, x (i")^ ■ The eigenvector, 
is a column vector with components CRL,i- Moreover, 1 is the unit matrix. 
El, is a diagonal matrix, and iJ is a Hermitian matrix. Vectors and diagonal 
matrices are denoted by lower-case Latin and Greek characters, and matrices 
by upper-case Latin characters. Exceptions to this rule are: Y (f ) , the vector 
of spherical harmonics, the site and angular-momentum indices (subscripts) 
R, L, I, and A, and the orders (superscripts) L, Af, and N. Operators are 
given in calligraphic, like Ti., and an omitted energy argument means that 
e = s^. 

With the (f) (r)'s being orthonormal in the ASA, the LMTO overlap matrix 
in (^) is seen to factorize to 1st order, and it is therefore simple to transform 
to a set of nearly orthonormal LMTOs: 



X(r) = x(i") 



1 



(7) 



{x\n-s,\x) 



H 



H 



1 + {H 



0(r) = 0(r) - 0(r)(0| 0) 
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Here, the energy-derivative function, (j) (r) , equals (r) , orthogonalized to 
(j) (r) . Finally, we may transform to a set of orthonormal LMTOs: 



X(r) = 



-1/2 



(8) 



H-e, 



We thus realize that of the Hamiltonians considered, H is of 1st, H is of 
2nd, and H is of 3rd order. As the order increases, and the energy window 
-inside which the eigenvalues of the Hamiltonian are useful as single-electron 
energies- widens, the real-space range of the Hamiltonian increases. For real- 
space calculations it is therefore important to be able to express 
a higher-order Hamiltonian as a power series in a lower-order Hamiltonian 
like in and (^), because such a series may be truncated when the energy 
window is sufficiently wide. 

The energy-derivative of the radial function ip (e, r) depends on the en- 
ergy derivative of its normalization. If we choose to normalize according to: 

(fi {e , r)'^ r'^ dr — 1, then it follows that (p (r) ip (r) r'^dr — 0. Choosing 
another energy-dependent normalization: ip{e,r) = (e, r) [1 + (e — e^) o] , 
specified by a constant o, then we see that: (p (r) = (f (r) + ip (r) o. Chang- 
ing the energy derivative of the normalization thus adds some Lp (r) to (p (r) 
and thereby changes the shape of the 'tail function' tp (r) . Since all LMTOs 
(0) should remain smooth upon this change, also H must change, and so 
must all LMTOs in the set. The diagonal matrix (^0 | 0^ , whose elements 

are the radial overlap integrals: o = <p> (r) ip (r) r^dr, thus determines the 
LMTO representation, and the first and the last equations (|7[) specify the 
linear transformation between representations. Values of the diagonal ma- 
trix (^(f) I 0^ exist, which yield short range for the Ist-order Hamiltonian H 

and, hence, for the LMTO set (|l|). Such an H is therefore a two-center TB 
Hamiltonian and such an LMTO set is a first-principles TB basis. 

In order to obtain an explicit expression for one needs to find the 
spherical-harmonics expansions about the various site for a set of smooth 
MTO envelope functions. For a MT-potential, which is flat in the intersti- 
tial, the envelope functions are wave-equation solutions with pure spherical- 
harmonics character near the sites. Consistent with the idea behind the ASA 
-to use 'space-filling spheres'- is the use of envelope functions with fixed en- 
ergy, specifically zero, which is a reasonable approximation for the kinetic 
energy between the atoms for a valence state. The envelope functions in the 
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ASA are thus screened multipole potentials, with the screening specified by 
a diagonal matrix of screening constants, am, related to the radial overlaps 
Ofl/. The expansion of a bare multipole potential at site R about a different 
site R' is well known: 



Here, Z" = Z' + Z and m" = m' — m. With suitable normalizations, the bare 
structure matrix, S°, can be made Hermitian. The screened structure matrix 
is now related to the bare one through a Dyson equation: 

(S")-^ = (S")"'- a, (9) 

which may be solved by inversion of the matrix S*^ — a^^ . This inversion may 
be performed in real space, that is in R- rather than in k-representation, 
provided that the screening constants take values known from experience to 
give a short-ranged S". 

In the end, it turns out that all ingredients to the LMTO Hamiltonian 

and overlap integrals, H, (^tp \ , and ^0 | t/)^ , may be obtained from the 

screened Korringa-Kohn-Rostoker (KKR) matrix in the ASA: 

^%L',B.L (e) ^ Pri (e) Sr'rSl'l - S%^,^j,j^. (10) 

Here, p*' (e) is a diagonal matrix of potential functions obtained from the 
radial logarithmic derivative functions, d{ip{e,s)} = c)ln (e, r) | /(9 lnr|^ , 
evaluated at the MT-radius, and p" (e) is related to p° (e) via the diagonal 
version of Equation (^. The results are: 



H = - K = £y-pp ^ + p ^Sp 2 = c + d2Sd2, 
= 2! =2!^' 3! ^3!i' 



(11) 



expressed in terms of the KKR matrix, renormalized to have K ^ \ : 

K{e) = k-iK{e)k-i = p (e) p^^ - p^^S p^i (12) 

This corresponds to the partial-wave normalization: ipir)^ r'^dr — 1, and 
K (e) is what in the 2nd-generation method is denoted —h (e) , but 

since the current notation identifies matrices by capitals, we have changed. 
The LMTO Hamiltonian and overlap matrices are thus expressed solely in 
terms of the structure matrix S and the potential functions p {e) , specifically 
the diagonal matrices p, p, p, and P . It may be realized that the nearly- 
orthonormal representation is generated if the diagonal screening matrix in 
(Bh is set to the value 7, which makes p''' vanish. 
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For calculations |^,^ 10 1 which employ the coherent-potential approxima- 
tion (CPA) to treat substitutional disorder, it is important to be able to 
perform screening transformations of the Green matrix: 

G"(z)^K"(z)-i = [p"(z)-S"]-\ (13) 

also called the resolvent, or the scattering path operator in multiple scattering 
theory ||ll|. In the 2nd generation MTO formalism, (e) was denoted g°' (e) . 
This screening transformation is: 

G^(z) = (/3-a)^ + ^G (14) 

and is seen to involve no matrix multiplications, but merely energy-dependent 
rescaling of matrix elements. As a transformation between the nearly or- 
thonormal, /3=7, and the short-ranged TB-representation, Eq. ( p^ has been 
useful also in Green- function calculations for extended defects, surfaces, and 
interfaces [^ P^l|l^ , [l3| , |l4| . However, calculations which start out from the 
unperturbed Green matrices most natural for the problem -namely those 
obtained from LMTO band-structure calculations in the nearly orthonor- 
mal representation for the bulk systems- have usually been limited to 2nd- 
order in z — because p''' (z) is linear to this order, and because attempts 
to use 3rd-order expressions for p'*' (z) employing the potential parameter 

p'' = 3! p'' (^(j) I 0^ 1 induced false poles in the Green matrix. 

What is not intelhgible in the TB-LMTO-ASA method is that the LMTO 
expansion (|l|) must include all L"s until convergence is reached throughout 
each sphere, and all i?"s until space is covered with spheres. This means 
that the LMTO-ASA basis is minimal -at most- for elemental, closely packed 
transition metals, the case for which it was in fact invented jl^. The supreme 
computational efficiency of the method soon made self-consistent density- 
functional calculations possible, and not only for elemental transition 
metals, but also for compounds. In order to treat open structures such as 
diamond, empty spheres were introduced as a device for describing the repul- 
sive potentials in the interstices ||l^. All of this then, led to misinterpretations 
of the wave-function related output of such calculations in terms of the com- 
ponents of the one-center expansions (|l|), typically the numbers of s, p, and 
d electrons on the various atoms (including in the empty spheres!) and the 
charge transfers between them. Absurd statements to the effect that CsCl is 
basically a neutral compound with the Cs electron having a bit of s-, more 
p-, quite some d-, and a bit of /-character were not uncommon. Many prac- 
titioners of the ASA method did not realize that the role of the MT-spheres 
is to describe the input potential, rather than the output wave- functions. For 
the latter, the one-center expansions truncated outside the spheres constitute 
merely a decomposition which is used in the code for selfconsistent calcula- 
tions. The strange Cs electron is therefore little more than the expansion 
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about the Cs site of the tails of the neighboring CI p electrons spilling into 
the Cs sphere. That latter MT-sphere must of course be chosen to have about 
the same size as that of CI, because only then is the shape of the Cs^Cl^ 
potential in the bi-partitioned structure well described. 

Now, the so-called high partial waves -they are those which are shaped 
like r' in the outer part of the sphere where the potential flattens out- do 
enter the LMTO expansion (0), but not the eigenvalue problem (||) or the 
equivalent KKR equation: 

K{e,)c,^0, (15) 

because they are part of the MTO envelope functions. This property of having 
the high-Z limit correct is a strength of the MTO method, not shared by for 
instance Gaussian orbitals, which are solutions of (^) for a parabolic potential. 
There are, however, also other partial waves -like the Cs s-waves, d-waves 
in non-transition metal atoms, /-waves in transition-metal atoms, s-waves 
in oxygen and fluorine, and in positive alkaline ions, and all partial waves 
in empty spheres- which for the problem at hand are judged to be inactive 
and should therefore not have corresponding LMTOs in the basis. In order 
to get rid of such inactive LMTOs, one must first -by means of (^) or (p^- 
transform to a representation in which the inactive partial waves appear only 
in the 'tails' (second term of (|l|)) of the remaining LMTOs; only thereafter, 
the inactive LMTOs can be deleted. This down-folding procedure works for 
the LMTO- ASA method, but it messes up the connection between the LMTO 
Hamiltonian (|^)-(|l^) and the KKR Green-function formalisms (|ll|)-([l5|), and 
it is not as efficient as one would have liked it to be E.g., the Si valence 
band cannot be described with an sp LMTO basis set derived by down-folding 
of the Si d- as well as all empty-sphere partial waves . 

The basic reason for these failures is that the ASA envelopes are chosen 
to be independent of energy -in order to avoid energy dependence of the 
structure matrix- because this is what forces us to carry out explicitly the 
integrals involving all partial waves in all spheres throughout space. What 
should be done is to include all inactive waves, ipi (e, r) , in energy- dependent 
MTO- envelopes, and then to linearize these MTOs to form LMTOs. This 
has been achieved with the development of the LMTO method of the 3rd- 
generation [ p^|j20[] , and will be dealt with in the present paper. The reason 
why energy linearization still works in a window of useful width, now that 
the energy dependence is kept throughout space, is due to the screening of 
the wave-equation solutions used as envelope functions [pit . 

As an extreme example, it was demonstrated in Fig. 7 of Ref. |^ -and 
we shall present further results in Fig. |l^ below- how with this method one 
may pick the orbital of one band, with a particular local symmetry and en- 
ergy range, out of a complex of overlapping bands. This goes beyond the 
construction of a Wannier function and has relevance for the treatment of 
correlated electrons in narrow bands P3,E3[. Another example to be treated 
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in the present paper is the valence and low-lying conduction-band structure of 
GaAs calculated with the minimal Ga spd As sp basis [|^. Other examples, 
not treated in this paper, concern the calculation of chemical indicators, such 
as the crystal-orbital-overlap-projected densities of states (COOPs) |^ for 
describing chemical pair bonding. These indicators were originally developed 
for the empirical Hiickel method where all parameters have been standard- 
ized. When one tries to take this over to an ab initio method, one immediately 
gets confronted with the problems of representation. For instance, COOPs 
will vanish in a basis of orthonormal orbitals. Therefore, the COOPs first had 
to be substituted by COHPs, which are Hamiltonian- rather than overlap pro- 
jections, but still, the LMTO-ASA method often gave strange results -for the 
above mentioned reasons [Q. What one has to do is -through downfolding- 
to chose the chemically- correct LMTO Hilbert space and -through screening- 
choose the chemically correct axes (orbitals) in this space. Only with such 
orbitals, does it make sense to compute indicators p7| , p8| . 

A current criterion for an electronic-structure method to be accurate 
and robust is that it can be used in ab initio density-functional molecular- 
dynamics (DF-MD) calculations |29). According to this criterion, hardly any 
existing LMTO method -and the LMTO-ASA least of all- is accurate and 
robust. 

Most LMTO calculations include non-ASA corrections to the Hamiltonian 
and overlap matrices, such as the combined correction for the neglected in- 
tegrals over the interstitial region and the neglected high partial waves. This 
brings in the first energy derivative of the structure matrix, S, in a way which 
makes the formalism clumsy . The code for the 2nd-generation LMTO 
method is useful ||3^ and quite accurate for calculating energy bands, be- 
cause it includes downfolding in addition to the combined correction, as well 
as an automatic way of dividing space into MT-spheres, but the underlying 
formalism is complicated. 

There certainly are LMTO methods sufficiently accurate to provide struc- 
tural energies and forces within density-functional theory [^2|j3^ , |3^ , |35| , ^6|j37|j38 
but their basis functions are defined with respect to MT-potcntials which do 
not overlap. As a consequence, in order to describe adequately the correspond- 
ingly large interstitial region, these LMTO sets must include extra degrees of 
freedom, such as LMTOs centered at interstitial sites and LMTOs with more 
than one radial quantum number. The latter include LMTOs with tails of dif- 
ferent kinetic energies (multiple kappa -sets) and LMTOs for semi-core states. 
Moreover, these methods usually do not employ short-ranged representations. 
Finally, since a non-overlapping MT potential is a poor approximation to the 
self-consistent potential, these methods are forced to include the matrix el- 
ements of the full potential. Existing full-potential methods are thus set up 
to provide final, numerical results at relatively low cost, but since they are 
complicated, they have sofar lacked the robustness needed for DF-MD, and 
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their formalisms provide little insight to the physics and chemistry of the 
problem. 

One of the early full-potential MTO methods did fold down extra or- 
bitals and furthermore contained a scheme by which the matrix elements of 
the full potential could be efficiently approximated by integrals in overlapping 



spheres |38 . The formalism however remained complicated, and the method 
apparently never took off. A decade later, it was shown |^,^ that the MT- 
potential, which defines the MTOs -and to which the Hamiltonian (Q) refers- 
may in fact have some overlap: If one solves the exact KKR equations 
with phase shifts calculated for MT-wells which overlap, then the resulting 
wave function is the one for the superposition of these MT-wells, plus an er- 
ror of 2nd order in the potential-overlap. This proof will be repeated in Eq. 
( p?!) of the present paper, and in Figs. |lj and |l^ we shall supplement the 



demonstration in Ref . 1 20 1 that this may be exploited to make the kind of ex- 
tra LMTOs mentioned above superfluous, provided that the MTO-envelopes 
have the proper energy dependence, that is, provided that 3rd generation 
LMTOs are used. Presently we can handle MT-potentials with up to ~60% 
radial overlap (s/? -I- sri < 1.6 |R — R'|), and it seems as if such potentials, 
with the MT-wells centered exclusively on the atoms, are sufficiently realistic 
that we only need the minimal LMTO set defined therefrom |20||4^ . It may 
even be that such fat MT-potentials, without full-potential corrections to the 
Hamiltonian matrix, will yield output charge densities which, when used in 
connection with the Hohenberg-Kohn variational principle for the total en- 
ergy will yield good structural energies |Q. Hence, we are getting rid of 
one of the major obstacles to LMTO DF-MD calculations, the empty spheres. 

Soon after the development of the TB-LMTO-ASA method, it was real- 
ized [ pl that the full charge density produced with this method -for cases 
where atomic and interstitial MT-spheres fill space well- is so accurate, that 
it should suffice for the calculation of total energies, provided that this charge 
density is used in connection with a variational principle. However, it took 
ten years before the first successful implementation was published ||45[] . The 
problem is as follows: The charge density, p (r) = J2T'^ I'^i > most simply 
obtained in the form of one-center expansions: 



(r) = X! X! / ImGflx.flX' (z) ipRL' [z, r^)* — , 



(16) 



where G{z) = K (x)"^ , as can be seen from (|l|) and @, but these expansions 
have terribly bad L-convergence in the region between the atoms and cannot 
even be used to plot the charge-density in that region. That was made possible 
by the transformation to a short-ranged representation, because one could 
now use: 



r,T r,i T / J OCC ^ 



RL R'L 



XR'L'{VR,)\ (17) 
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where the L-sums only run over active values, and where the double-sum over 
sites converges fast. Nevertheless, to compute a value of xb.l (r) with r far 
away from a site, one must evaluate the LMTO envelope function, which is a 
superposition of the bare ones, Yl {^r) lf''R^ ^ and this means that ( p^ ) actu- 
ally contains a 4-double summation over sites. At that time, this appeared to 
make the evaluation of p (r) at a sufficient number of interstitial points too 
time-consuming for DF-MD, although the full charge density from (jl^) was 
used routinely for plotting the charge-density, the electron-localization func- 
tion a.s.o. In order to evaluate the total energy, the full charge density 
must also be expressed in a form practical for solving the Poisson equation. If 
one insists on a real-space method, then fast Fourier transformation is not an 



option. In Fig. 12 of the present paper, we shall present results of a real-space 



scheme |47| , t48|| used in connection with 3rd-generation LMTOs for the phase 
diagram of Si |49|. This scheme is presently not a full-potential, but a full 
charge-density scheme, and the calculation of inter-atomic forces has still not 
been implemented. 

With 3rd generation LMTOs |ll9|,|o|, the simple ASA expressions (0)- 
still hold, provided that (e, r) is suitably redefined, and that K (e) 
is substituted by the proper screened KKR matrix whose structure matrix 
depends on energy. The LMTO Hamiltonian and overlap matrices are given 
in terms of and its first three energy derivatives, K, K, and K, which 
are not diagonal. Downfolding, the interstitial region, and potential-overlap 
to first order are now all included in this simple AS A- like formalism Q. In 
due course, we thus hope to be able to perform DF-MD calculations with an 
electronic Hamiltonian which is little more complicated than (||), (Q), or (||). 

A final problem with the LMTO basis is that even with the conventional 
spd-basis and space-filling spheres, the LMTO set is insufficient for cases 
where semi-core states and excited states must be described by one minimal 
basis set, and in one energy panel. This problem becomes even more acute 
in the 3rd-generation method where, due to the proper treatment of the 
interstitial region, the expansion energy must be global, that is, is 
now the unit matrix times e^, rather than a diagonal matrix with elements 
SuRiSbb/Slli . The same problem was met when attempting to apply the 
formally elegant relativistic, spin-polarized LMTO method of Ref. ||5^ to 
narrow, spin-orbit split /-bands. Finally, as MT-spheres get larger, and as 
more partial waves are being folded into the MTO envelopes, the energy 
window inside which the LMTO basis gives accurate results shrinks. This 
means, that the 3rd-generation LMTO method described in |Q may not be 
sufRciently robust. 

The idea emerging from the LMTO construction (|l|) seems to be: Divide 
space into local regions inside which Schrodinger's equation separates due to 
spherical symmetry and which are so small that the energy dependence of 
the radial functions is weak over the energy range of interest. Then expand 
this energy dependence in a Taylor series to first order around the energy e^, 
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at the center of interest: (f) (e, r) « (r) + (e — Eu) 4> (r) • Finally, substitute 
the energy by a Hamiltonian to obtain the energy- mdependent LMTO. The 
question therefore arises (Fig. ^): Can we develop a more general, polynomial 
MTO scheme of degree TV, which allows us to use an iVth-order Taylor series 
or -more generally- allows us to use a mesh of iV + 1 discrete energy points, 
and thereby obtain good results over a wider energy range, without increasing 
the size of the basis setl Such an NMTO scheme has recently been developed 
[ pT| and shown to be very powerful |Q. We shall preview it in the present 
paper. 





Fig. 1. Quadratic approximation to the energy dependence of a partial wave for a 
condensed (Taylor) and a discrete (Lagrange) mesh. 



Most aspects of the 3rd-generation LMTO method have been dealt with in 
a set of lecture notes Q and a recent review ||20j. Here, we shall try to avoid 
repetition but, nevertheless, give a self-contained description of two selected 
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aspects of the new method: the basic concepts and the new polynomial NMTO 
scheme, to be presented here for the first time. 

We first explain (Sect. |^) what the functions (/)(e, r) actually are in the 
3rd-generation formalism. This we do using conventional notation in terms of 
spherical Bessel functions and phase shifts -like in Ref . pl| |- and only later, 
we renormalize to the notation used in Refs. and ]20|. It turns out that 
the bare 0's are the energy- dependent MTOs of the 1st generation |52|. The 
screened (/)'s are the screened, energy-dependent MTOs of the 2nd generation 
[ pT| , with the proviso that = e. This proviso -together with truncations of 
the screening divergencies at the sites, inside the so-called screening spheres- 
is what makes the screened 0's equal to the so-called unitary |ig[| or kinked 
partial waves in the formalism of the 3rd generation. We then derive 
the screened KKR equations and repeat the proof from Refs. ||2l[ and 
that overlapping MT-potentials are treated correctly to leading (1st) order in 
the potential overlap. Towards the end of this first section, we introduce the 
so-called contracted Green function (e, r) G (e) , which will play a crucial 
role in the development of the polynomial NMTO scheme, and we derive the 
3rd-generation version of the scaling relation ( [T^ ) for screening the Green 
function. 

In Sect. I we show how to get rid of the energy dependence of the kinked- 
partial wave set: First, we introduce a set of energy-dependent NMTOs, 
X^^^ {£, r) , which -hke the </) (e, r) set- spans the solutions of Schrodinger's 
equation for the chosen MT-potential, and whose contracted Green function, 
X*-^-* {e,r)G{e) , differs from (e, r) G (e) by a function which is analytical 
in energy. Like in classical polynomial approximations, we choose a mesh of 
arbitrarily spaced energies, eo,...,eAr, and subsequently adjust the analyti- 
cal function in such a way that, x^^"* (^o^r) = ... = x^^"* (^JV,r) • The latter 
then, constitutes the set of energy- mdependent NMTOs. The Oth-order set, 
T^(o) (i-) ^ ig seen to be the set of kinked partial waves, (j) {eq, r) , at the energy 
Eo, and the Ist-order set, x^^'' (r) , to be the set of tangent or chord-LMTOs 
-depending on whether the mesh is condensed or discrete. For the case of a 
condensed mesh -which is the simplest- the matrices, which substitute for 
the energies in the Taylor series (|l|) -generalized to A^th order- turn out to 



be: 




for 1<M< N, 



(18) 
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in terms of the Mth and the {M — l)st energy derivatives of the Green matrix. 
Moreover, the expressions for the Hamiltonian and overlap matrices are: 



(AT) 

iV! 



(27V) 

G 
(2iV)! 



m 



(19) 



G 



Nl 



(2Ar+l) 

G 

(2iV+ 1)! 




which, for = 1, are easily seen to reduce to (||) upon insertion of ([ll]). In 
retrospect, it is convenient that these basic NMTO results are expressed in 
terms of energy derivatives of the Green matrix G (e) -rather than in terms of 
those of its inverse ,the KKR matrix, as we are used to from the LMTO-ASA 
method (|rT|)~ because if we imagine generalizing (^ to A^th order and using 
it to form the Hamiltonian and overlap matrices like in ^ , then each matrix 
will consist of N'^ terms, among which a number of relations can be shown to 
exist. We also realize, that the problem mentioned above about using Green 
matrices beyond 2nd order in z — is solved by using -instead of G (z)- the 
NMTO Green function: 



AN) 



(JV) 

(W)\"^ ^ _^ 

/ iV! 



{2N) 

G 
(27V)! 



- (z - Si,) 



(2N+1) 

G 

(27V + 1) 



(JV) 



which equals G (z) to {2N + l)st order. This Green function has the addi- 
tional advantage of allowing for a simple treatment of non-MT perturbations. 
We admit that this route to energy-independent MTO basis sets has little 
in common with the twisted path we cut the first time, but once found, it is 
easy to accept and understand the results -which are simple. 

In practice, it is cumbersome to differentiate a KKR matrix -not to speak 
of a Green matrix- many times with respect to energy. Hence, one uses a dis- 
crete energy mesh. With that, the derivatives in (|l8| ) and the pre- and post 
factors in ( |l9| ) and (^0|) turn out to be divided differences, while those at 
the centers of (^9|) turn out to be the highest derivative of that approximat- 
ing polynomial which is fitted not only to the values of G (e) at the mesh 
points, but also to its slopes. Hence, they are related to classical Hermite 
interpolation p3| . 

In both Sections || and ^, special attention is paid to the so-called triple- 
valuedness, because this was not previously explained in any detail, but has 
turned out to be crucial for the further developments and will be even more 
so when we come to evaluate the inter-atomic forces. A related aspect is the 
fact that a screening transformation in the formalism of the 3rd-generation is 
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linear as regards the envelope functions, but non-linear as regards the NM- 
TOs. This means, that changing the screening, changes the NMTO Hilbert 
space. This was not the case for 2nd-generation LMTOs. This is the reason 
why we took care to denote the nearly-orthonormal and orthonormal LMTO 
sets arrived at by the linear transformations (Q) and by respectively x 
X, rather than by x'' and x'^i as in the 2nd-generation LMTO scheme, where 
screening transformations were linear and denoted by superscripts. Screening 
transformations like (^) and ( jl^ ) still hold for the 3rd-generation structure- 
and Green-matrices, but the partial waves providing the spatial factors of 
the Green function (seejl^)) are different: they have tails extending into the 
interstitial region. A tail is attached continuously, but with a kink, at the 
screening sphere, which is concentric with, but smaller than, its own MT- 
sphere, and the resulting kinked partial wave, or Oth-order energy-dependent 
MTO, is -for the purpose of evaluating its properties in a simple, approxi- 
mate way- triple-valued in the shell between these two spheres. The radii, 
o-RL, define the screening and determine the shape of the MTO envelopes. 
Now, for a superposition of kinked partial waves given by a solution of the 
KKR equations the kinks and the triple- valuedness cancel, but for a 

single NMTO, a triple-valuedness of order (r — a)^^^^ (e^ — Eq) ... (e^ — En) 
-which is the same as the error caused by the energy interpolation- remains. 
For this reason: The smaller the screening radii -i.e. the weaker the screening- 
the smaller the energy window inside which an energy- mdependent NMTO 
set gives good results. The extreme case is the bare (a 0) N — set, 
which is the set of Ist-generation MTOs [Q, but defined without freezing 
the energy dependence outside the central MT-sphere. The tail-cancellation 
condition for this set leads to the original KKR equations which -we 
know- must be solved energy-by-energy, that is, the energy window can be 
very narrow, depending on the application. Specifically, for free electrons the 
width is zero. 

At the end of Sect. ^, we demonstrate the power of the new NMTO 
methods by applying the differential and discrete LMTO, QMTO, and CMTO 
variational methods to the valence and conduction-band structure of GaAs 
using a minimal Ga spd As sp basis, and to the conduction band of CaCu02 
using only one orbital, all others being removed by massive downfolding 
[ p4| . We also give simple expressions for the charge density and show the 
total energy as a function of volume for the various crystalline phases of Si 
calculated with the full-charge, diflFerential LMTO method |4^,||,|9|. Finally, 
numerical results are presented for the error of the valence-band energy of 
diamond-structured Si -as a function of the potential overlap- obtained from 
LMTOs constructed for a potential whose MT-wells are centered exclusively 
on the atoms. In addition, results of a scheme which corrects for the error of 
2nd order in the overlap will be presented |4^ . 

In Sect. ^ we show that energy-dependent, linear transformations of the 
set of kinked partial waves -such as a normalization- merely leads to similar- 
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ity transformations among the NMTO basis functions and, hence, does not 
change the Hilbert space spanned by the NMTO set. 

This is exploited in Sect. || to generate nearly orthonormal basis sets, 
X^^'^ (r) , for which the energy matrices defined in ( p^ ) become Hermitian, 
Hamiltonian matrices, Z/*^*^^. We also show how to generate orthonormal 
sets, x''^'' (i") J of general order, and we demonstrate by the example of the 
minimal MTO set for GaAs that this technique works numerically efficiently 
-at least up to and including = 3. This development of orthonormal basis 
sets should be important e.g. for the construction of correlated, multi-orbital 
Hamiltonians for real materials [^3y54[| . 

In the last Sect. ^ we show explicitly how -for iV = 1 and a condensed 
mesh- the general, nearly-orthonormal NMTO formalism reduces to the sim- 
ple ASA formalism of the present Overview. 

In the Appendix we have derived those parts of the classical formalism for 
polynomial approximation -Lagrange, Newton, and Hermite interpolation- 
needed for the development of the NMTO method for discrete meshes . 

2 Kinked partial waves 

In this section we shall define Oth-order energy-dependent MTOs and show 
that linear combinations can be formed which solve Schrodinger's equation for 
the MT-potential used to construct the MTOs. The coefficients of these linear 
combinations are the solutions of the (screened) KKR equations. By renor- 
malization and truncation of the irregular parts of the screened MTOs inside 
appropriately defined screening spheres, these Oth-order energy-dependent 
MTOs become the kinked partial waves of the 3rd generation. 

If we continue the regular solution ipi^i (e, r) of the radial Schrodinger 
equation (^) for the single potential well, vr (r) , smoothly outside that well, 
it becomes: 

(pm{e,r) = ni{Kr) - ji{Kr) cot r]Ri{e) = ip°jii{e,r) , for r > sr, (21) 

in terms of the spherical Bessel and Neumann functions, ji {nr) and ni (nr) , 
which are regular respectively at the origin and at infinity, and a phase shift 
defined by: 

^ n (ks) dln\ip {e,r)\ /dliir\^ — dhi\n {Kr)\ /dhir\^ 

^ j {ns) dln\(p {s,r) \ /d\nr\^ — 9 In \j {Kr) \ /91nr|^ ' 

In the latter expression, we have dropped the subscripts. Note that we no 
longer distinguish between 'inside' and 'outside' kinetic energies, e — v{r) and 
= s — Vmtz , and that we have returned to the common practice of setting 
Vmtz = 0. If the energy is negative, ni [nr) denotes a spherical, exponentially 
decreasing Hankel function. Note also that -unlike in the ASA- the radial 
function is not truncated outside its MT-sphere, and is not normalized to 
unity inside. In fact, we shall meet three different normalizations throughout 
the bulk of this paper, and ( |2l| ) is the first. 
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Fig. 2. Bare Si p MTO according to Eq.(||) 



2.1 Bare MTOs 

The bare, energy-dependent muffin-tin orbital (MTO) remains the one of the 
1st generation fs^ : 

(jiRL (e, r) = Yl (?) [ipm (e, r) + ji {nr) cot i^ri (e)] 

^ In; (rer) for r > sr 

= Yl (r) [ipRi (e, r) - <^?j, (e, r) («r)] , (22) 

and is seen to have pure angular momentum and to be regular in all space. The 
reason for denoting this Oth-order MTO 0(e,r) , rather than ;i^(^=o) (e, r) , 
should become clear later. 

In Fig. ^ we show the radial part of this MTO for a Si p-orbital, a MT- 
sphere which is so large that it reaches 3/4 the distance to the next site in 
the diamond lattice, and an energy in the valence-band, which -in this case 
of a large MT-sphere- is slightly negative (see Fig. 11 in Rcf. |Q). The full 
line shows the MTO as defined in (p2|), while the various broken lines show 
it 'the 3-fold way': The radial Schrodinger equation for the potential v (r) is 
integrated outwards, from the origin to the MT radius, s, yielding the regular 
solution, (p(e,r), shown by the dot-dashed curve. At s, the integration is 
continued with reversed direction and with the potential substituted by the 
flat potential, whose value is defined as the zero of energy. This inwards 
integration results in the radial function 'seen from the outside of the atom', 
ip° (e, r) , shown by the dotted curve. The inwards integration is continued to 
the origin, where ip° (e, r) joins the 'outgoing' solution for the flat potential, 
that is the one which is regular at inflnity: n {nr) . The latter is the envelope 
function for the bare MTO. 
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As usual, the envelope-function for the MTO centered at R may be ex- 
panded in spherical-harmonics about another site R' (7^ R): 

Km {ktr) Yl (fi?) ^^ji' {KrR,)YL' {rR')BR,L',RL (e) , 

L' 

where the expansion coefficients form the Hermitian KKR structure matrix: 
BwL',RL (e) = ^47rr'+''-'"CLL.z.- Krir {k |R- R'|) Y;,^^„ (r^') 

(23) 

as conventionally |^ defined, albeit in i?-space. The spherical harmonics are 
as defined by Condon and Shortley, m" = m! — m, the summation runs over 
I" = \V - l\ , \V - ?| + 2, V -I- /, and is real, because Cll'L" = 

^Yi^{f)Yl.(r)Yi^„{f)df. 

If for the on-site elements of B (e) , we define: Brl rl' (e) = 0, and use the 
notation: (e, tr) = /; (ktr) Yl (tr) , as well as the vector-matrix notation 
introduced in connection with we may express the spherical-harmonics 
expansion of the bare envelope about any site as: 

KTi (e, r) = j (e, r) B (e) -1- nn (e, r) . (24) 

When we now form a linear combination, ^R^fj^RL {£,''^r) crl, of energy- 
dependent MTOs (p^, and require that it be a solution of Schrodinger's 
equation, then the condition is that, inside any MT-sphere (i?') and for any 
angular momentum {L') , the contributions from the tails should cancel the 
ji' (nr) cot -qR' I' (e)-term from their own MTO, (jiR'u {£, Ti?/), thus leaving be- 
hind the term fR'i' (e, r) , which is a solution by construction. This gives rise 
to the original KKR equations : 

[Br'L'.RL (Si) + KCOt r]Ri {Si) SrirSl'l] CRL,i 

RL 

= ^ Kr'L',rl (ei) CRL,i = 0, (25) 

RL 

which have non-zero solutions, CRL,i, for those energies, £i, where the deter- 
minant of the KKR matrix vanishes. 

With those equations satisfied, the wave function is 

00 I' 

'^(f'RL {£i,rR)cRL,i = X! X! (jR')cRfL'S + (26) 

RL l'=Om' = -l' 

t^-R' ~ '^Rl ^«)] (f fl) CRL,^ 

R^R' L 

near site R'. Since according to ( pi] ) the function ~ (p° vanishes outside its 
own MT-sphere, the terms in the second line vanish for a non-overlapping 
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MT-potential so that, in this case, ( [2q ) solves Schrodinger's equation exactly. 
If the potential from a neighboring site (i?) overlaps the central site {R'), then 
fRL — f°RL tongues stick into the MT-sphere at R' . The radial part of such a 
tongue is \ {sr - r^f vr {sr) (Prl {sr) , to lowest order in sr - tr, as may 
be seen from the radial Schrodinger equation (j^). Let us now operate on 
the smooth function S'i (r) = J2rl 'I^RL (sii'^R) CRL,i , of which ( p^ ) is the 
expansion around site i?', with Ti. — Si as given by (^) to find the error: 

VR' (rfl/) Y ^R) ~ ^°Rl i^^lTR)] (ffl) CRL^^ (27) 

R' R^R' L 

^ pairs 
RR' 

This shows that the wave function ( ^61 ) solves Schrodinger's equation for the 
superposition of MT-wells to within an error, which is of second order in the 
potential overlap ||2l| , p0t . 



{sR' - rR, f + {sR - rRf VR (sr) (r) . 



2.2 Screened MTOs 

Screening is the characteristic of 2nd-generation MTOs and was first dis- 
covered as the transformation (^ to a ne arly - orthonormal representation, in 
which the Hamiltonian is of second order fSq , ^ . Shortly thereafter it was re- 
alized that there exists a whole set of screening transformations which may be 
used to make the orbitals short ranged^ so that the structure matrix may be 
generated in real space. It was also realized that the screening transformation 
could be used to downfold inactive channels and, hence, to produce minimal 
basis sets These applications were all for the ASA with k^=0. Only 

long time after ||2l|| , did it become clear that screening would work for positive 
energies as well, and at that time a fourth virtue of screening became clear, 
namely, that sceening the range of the orbitals, simultaneously reduces their 
energy dependence to the extent that the full energy dependence may be kept 
in the interstitial region, thus making the K^=0-part of the ASA superfluous. 
Most of this was shown in the last paper on the 2nd-generation formalism 
[ pT| . Nevertheless, this paper was unable to devise a generally useful recipe 
for choosing the energy-dependent screening constants, it failed to realize 
that screening allows the return to: n^—e, and for those reasons it missed the 
elegant energy-linearization of the MTOs achieved by the 3rd generation. 

The screened envelopes of the 2nd-generation method are linear superpo- 
sitions. 



(e, v) = n (e, r) 5" (e) , 



(28) 
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of the envelope functions, n(e,r), with the property that the spherical- 
harmonics expansions of the set of screened envelopes be: 

nn (e, r) S*" (e) = Kn" (e, r) = j° (e, r) (e) + ku (e, r) , (29) 

which are ( p^ with the substitutions: 

ji{Kr) -> (e,r) = j; (Kr) - n; (Kr)tanQ!fl,/„ (e) , (30) 

and: B (e) ^ B" (e) , which will be determined below. In contrast to its bare 
counterpart, a screened envelope does not have pure angular momentum, i.e., 
cannot be factorized as a radial function times a spherical harmonics, and 
it depends explicitly on its surroundings. The background phase shifts a (e) 
-which may even depend on m (see for instance Fig. |ll|)- specify the shapes 
of the screened envelopes. Whereas the bare envelopes are regular in all space 
-except at their own site where they diverge like F;™ (f) /r'+^- the screened 
envelopes diverge at any site where there is a finite background phase shift 
in at least one L-channel. 

Note that only in the Overview did we use ASA K^=0-notation with Greek 
letters denoting screening constants and S" the structure matrix. In the bulk 
of the present paper, we use Greek letters to denote background phase shifts, 
and B°' and 5*" to denote respectively the structure matrix and the screening 
transformation. 

We now find the screened structure matrix and the transformation matrix 
by expanding also the bare envelope on the left hand side of ( p9| ) by means 
of (|2^). Comparisons of the coefficients to /tn^/ (e, yri) and jl' (e, ffl') yield 
respectively: 

^° (e) = i_ ^^'^"(^) g" (g)^ and: (e) = B (e) 5" (e) (31) 

with the quantities regarded as matrices, e.g. tana is considered a diag- 
onal matrix with elements tsxiauL ^rr'^lu ■ As a result of (^l]): 

B-{e)-'^B{e)-' + '-^^, (32) 

which shows that, like the bare structure matrix, also the screened one is 
Hermitian. In contrast to the bare structure matrix, the screened one has 
non- vanishing on-site elements. For background phase shifts known to give a 
short-ranged i?" (e) , the inversion of the matrix B {e) n cot a (e) , implied 
by (|3^), may be performed in real space, although the hare structure matrix 
is long-ranged. Eq. (p2h is the k^=£ equivalent of the ASA 'Dyson equation' 



For the inactive channels {RL = /) , we choose the background phase 
shifts to be equal to the real phase shifts: 



ai (e) = rji (e) 



(33) 
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SO that for these channels, 

jfie,r) = ji (kv) - ni {Kr)tanrji (e) = -<^j (e, r) tan?]/ (e) . 

That is, we shape the set of screened envelope functions in such a way 
that, for the inactive channels, the radial functions, <pj , may be substi- 
tuted smoothly by the regular solutions, ipi (e, r) , of the radial Schrddinger 
equation. This is what we call downfolding. This substitution makes the 
screened envelopes become the so-called screened spherical waves, ip, of the 
3rd-generation method. Only the screened spherical waves corresponding to 
the remaining, so-called active channels (RL = A) will be used to construct 
the MTO; they are: 

^P^^ie^rn) ^ n%L{e,rn) + (34) 

J2 (e, rn') - {e, rj,,)] ^.^^^Yj {rn') B^^^ (e) , 

^ — ' K 

I 

which -in contrast to (e, rji)- are regular in all inactive channels, albeit 
irregular in the active channels. In (^), / = R'U. Below, we shall choose 
to truncate the active channels inside their screening spheres. Due to the 
augmentation (substitution) , the screened spherical waves do not transform 
linearly like (P8|). 

For the partial waves of high I, the phase shifts vanish due to the domi- 
nance of the centrifugal term over the potential term in the radial Schrodinger 
equation (^) . As a consequence, the matrices involved in the Dyson equation 
( p2[ ) -whose indices run over all active as well as inactive channels- truncate 
above a certain I of about 3-4. 

Before specifying our choice of background phase shifts for the active 
channels, let us define the energy-dependent, screened MTO analogous to 
the third equation ( p2[ ) as the (augmented) envelope function, plus a term 
proportional to the function 93 — 93°, which vanishes (quadratically) outside 
the central MT-sphere and has pure angular-momentum character. That is: 

(^RL (£> ^r) = Yl (rR) [(pRi (e, rn) 

and RL E A. Here, the coefficient to ip — (p° has been chosen in such a 
way that, in its own channel and outside any other MT-sphere, the screened 
MTO is ip" -I- j" cot rj" plus a term from the diagonal element of the screened 
structure matrix. 

To check this, we project onto the 'eigen-channel,' making use of (|4|), 
(p9|), (pi]), and (|30|), and neglecting any contribution from tpi (e, r/j')'s from 



-^^^^'^'^^^ tan,j^,(s) +^^^("'^^) 
-pins,rn)]+rRLie,rR) (35) 
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overlapping neighboring MT-spheres: 

r / r» \ 1 tan 77 

= Yp — n + [] + n tan a) cot 7;J h 77 + j 

tan 77" K 

„ „ „ tan 77 — tan a „ B" 

tan 77" K 

= + cot r/" + f — (36) 

For simplicity, we have dropped all arguments and indices in the last three 
lines. We see that the new phase shift, 77", is given by: 

tan77gx (e) = tan7yij; (e) - tana^x (e) , (37) 

as expected for the phase shift on the background of a. This is the same 
transformation as the one obtained from (H) for -B" [ey^ . The definition 
of the renormalized free radial solution given in (|35| ) may be written as: 

(e, r) ^ m {kt) - {s, r) cot ^rl (e) (38) 
= [?7,i [kt) i'AnriRi (e) - ji {nr)] cot rj'^]^ (e) , 

and ip'^i {e,rR) is the solution of the radial Schrodinger equation, normalized 
in such a way that it matches onto ip]^1 (e, r) at the MT radius, sr. The 
definition ( p8| ) reduces to ( pi] ) when a — 0. 

The set of screened MTOs now consists of the screened MTOs ( |35| ) of 
all active channels. Since the ip — ip° function has pure angular-momentum 
character, the mixed character of the screened MTO stems solely from the ip- 
function. The result of projecting the screened MTO onto an active channel 
R'L' different from its own is seen from ( p9[ ) to be: 

Vr^l'^rl (e,rij) = VR^L'rRL {s,rR)^r^,^, {e,rR,) ^^'^''^^ (39) 

when vri is so small that r lies inside only one MT-sphere, the one centered 
at R' . From ( ^ ) and ( |3^ ) it is then obvious that, in order to get a smooth 
linear combination 0^ (e, r^) of screened MTOs, all j"^ -functions must 
cancel. This leads to the condition that the energy must be such that the 
coefficients can satisfy 

E ^^A'A (e.) + cot Sa'a] cl, ^ J2 ^A'A cS,. = 0, (40) 

A A 

for all active R'L' = A' . These are the screened KKR equations, and if" (e) 
is the screened KKR matrix. If these equations are satisfied, the linear com- 
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bination of screened MTOs is: 

oo 

Y,rA {e^,VR) cX, ^Y.Y. VR'L' {e,,rR,)YL> (fflO c^'L',. + (41) 

A V=Q m' 

R^R' L 

near site i?'. As long as the MT-spheres do not overlap, this is a solution of 
Schrodinger's equation for the MT-potential and, if the potentials overlap, 
then the Lp — (f° tongues from the neighboring sites in the second line of 
(^ij) make the wave function correct to first order in the overlap |20| . This is 
exactly as in (|2^). The summation over spherical- harmonics around the cen- 
tral site includes the contributions —^pi {e, rff) tanry/ (e) J2a^? A (g) j 
provided by the screened-spherical-wave part of the MTO (see ( psj ) and (p^)). 

Although energy-dependent MTO sets with different screenings are not 
linearly related, they all solve Schrodinger's equation for the MT-potential 
used for their construction via the corresponding KKR equation. E.g. had one 
chosen a representation in which a channel making a significant contribution 
to a wave function Wi (r) with energy Ei = e were downfolded, then the cor- 
responding solution of the KKR equation (|4^) would arise from B" (e) being 
long ranged and, as a function of e, going through a zero-pole pair near e^. If 
the energy were now fixed at some energy e^, and the energy- independent set 
0" (£,y, r) were used as the Oth-order MTO basis in a variational calculation, 
then a useful result could in principle be obtained, but only if Si, were chosen 
very close to e^. 

Si p kinked partial wave 



0.4 



0.2 



0.0 

Si <-a a^s Si <-a 

Fig. 3. Si Pill member of a screened spd-set of Oth-order MTOs (see text and 
Eqs.(||,(g|-(||)). 
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2.3 Hard-sphere interpretation and redefinitions 

We now wish to choose the background phase shifts for the active channels in 
a way which reduces the spatial range and the energy dependence of the MTO 
envelopes. It is obvious, that for the orbitals to be localized, they must have 
energies below the bottom of the continuum of the background -defined as the 
system which has the same structure as the real system, but has all phase 
shifts equal to those of the background. Hence, the active a (e)'s should be 
defined in such a way that the energy band defined by: (e) + k cot a{e)\ = 
0, he as high as possible. 

The discovery of a useful way of determining this background, turned 
out to be the unplanned birth of the 3rd MTO generation Realiz- 
ing that the weakest point of the ASA was its solution of Poisson's -and 
not Schrodinger's- equation, and unhappy with the complexities of existing 
full-potential schemes, we ]57[ | were looking for those linear combinations of 
Hankel functions -like (^^- which would fit the charge density continuously 
at spheres. With Methfessel's formulation What we wanted was those 
solutions of the wave equation which are (jr) at their own sphere and for 
their own angular momentum, and zero at all other spheres and for all other 
angular momenta. This set was therefore named unitary spherical waves. The 
solution to this boundary-value problem is of course a particular screening 
transformation (^2|). 

Our way of defining the background was thus in terms of hard screening- 
spheres for the active channels; the larger the screening spheres, the larger the 
excluded volume and the higher the bottom of the continuum. The screening 
spheres are not allowed to overlap -at least not if all Z-channels were active, 
because then a unitary spherical wave would be asked to take both values, 1 
and 0, on the circle common to the central and an overlapping sphere. As a 
consequence, in order to reduce the range and the energy dependence of the 
MTO envelope functions, the screening spheres should in general be nearly 
touching. Now, since the screening radii, , control the shapes of the envelopes, 
the relative sizes of the screening spheres should be determined by chemical 
considerations, i.e. the a's may be covalent- or ionic radii in order that results 
obtained from an electronic-structure calculation be interpretable in terms of 
covalency, ionicity etc. Referring to the discussion in the Overview, one could 
say: The MT-spheres (s) are potential-spheres and the screening-spheres (a) 
are charge-spheres. 

Inspired by Ref . , practitioners of multiple-scattering theory -who tra- 
ditionally take the Kohn-Rostoker Green-function point of view- found 
another useful way of determining the background phase shifts, namely in 
terms of repulsive potentials |5^ . 

For a given active channel {RL = A), the radial positions, r — oa (e) , of 
the nodes of the background functions j" given by (^0[) are the solutions of 
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the equation: 

Ja (e, aA (e)) = ji {kua (e)) - rn {kqa (s)) tana^ (e) . 

Whereas attractive potentials usuaUy do not give positive radii -for an ex- 
ample, see the dotted curve in Fig. ^ repulsive potentials do, as may be 
seen from the radial Schrodinger equation (|^). For a hard repulsive potential, 
the position of the node is independent of energy and of I. What we shall 
use for the active channels are therefore screening-sphere radii, aA, which are 
independent of energy and which usually depend little on L among the active 
channels. In terms of such a screening radius, the corresponding background 
phase shift is given by: 

tana^ (e) = ji (naA) /ni (naA) ■ (42) 



Now, instead of having screened spherical waves (^4|) and MTOs ( p5[) 
whose active channels are irregular at the origin -the irregularities of the 
mactive channels were already gotten rid of by downfolding, followed by 
<Pj(£,r) —^ ipi{e,r) substitutions- we prefer that the active channels have 
merely kinks. This is achieved by truncating all active j"-functions inside 
their screening spheres, that is, we perform the substitution: 

, l""" (43) 

yji(Kr) - ni(Kr)ji{K.aA) Ini^KaA) for r > aA 

which is continuous but not diffcrentiable, for the screened spherical waves 
and for its own j"-function of the MTO -that is the second term on the 



last two lines of (36). With that substitution, a screened spherical wave, 
V'Si (^' ^r) ' vanishes inside all screening spheres of the active channels - 
except inside its own, where it equals ni {nrji) {^r) ■ This may be seen 
from ( |39| ) and the two first lines of (^) . Finally, if we renormalize according 
to: 

^RL (£> ^r) = ^RL (e, i-fl) /ni (naRL) (44) 

-note the difference between the superscripts a and a- we finally arrive at 
the screened (unitary) spherical wave as defined in Refs. [|l9[|20|1 . 

■0^^ {e, rn) is that solution of the wave equation which is Yl (f _r) on 
its own screening sphere, has vanishing Yl' (fi?')-average on the screening 
spheres of the other active channels, and joins smoothly onto the regular 
solutions of the radial Schrodinger equations of the inactive channels. In 
those, the regular Schrodinger solutions are, in fact, substituted for the wave- 
equation solutions. 

It is now obvious, that overlap of screening spheres will cause complicated, 
and hence long-ranged spatial behavior of the screened spherical waves, and 
the worse, the more spherical harmonics are active. 
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With the normahzation (p^, there is apparently no need for functions, 
like spherical Bessel and Neumann or Hankel functions, which have a branch- 
cut at zero energy, and this was the point of view taken in the first accounts 
[ p^[]20| of the 3rd-generation method. However, the normalization (Q) is not 
appropriate for a—0, and expressing the screened structure matrix in terms of 
the bare one (|2^) -which is the only one computable in terms of elementary 
functions- was slightly painful in Ref. p9| ]; moreover, in that paper down- 
folding was not presented in its full generality. In these respects, the present, 
conventional derivation is simpler, but it takes more equations. 

With the a ^ a redefinitions (^)-(|4^), the MTO remains as defined by 
35|) , but with the screened spherical waves and its own j"-fimction truncated 
as described above. We may also renormalize the MTO like in (^^: 

(t>RL (£> ^r) = (t>RL (£> ^r) /ni {naiiL) , (45) 

whereby these energy-dependent Oth-order MTOs become identical with the 
kinked partial waves of Refs. [|l9|,p0|. This normalization corresponds to: 



^?j?(e,afli)^l. (46) 

Note that this will cause the normalization of the radial Schrodinger-equation 
solution, if"" (e, r) , to depend on m in case the corresponding screening radius 
is chosen to do so. 

In Fig. H we show the screened counterpart of the bare Si p orbital in Fig. 
Ij. Since only the two first terms of (^5|) -but not the screened spherical wave- 
has pure angular momentum, we cannot plot just the radial wave function like 
in Fig. ||. Rather, we show the MTO together with its three parts along the 
[lll]-line between the central atom and one of its four nearest neighbors in 
the diamond structure. The positions of the central and the nearest-neighbor 
atoms are indicated on the axis (Si) , and so is the intersection with the cen- 
tral MT-sphere (s). The p orbital chosen is the one pointing along this [111] 
direction. The Si spd channels were taken as active, and to have one and 
the same screening radius, a — 0.75t, where t is half the nearest-neighbor 
distance, i.e., the touching-sphere radius. The places where the central and 
the nearest-neighbor screening spheres intersect the [lll]-line are indicated 
by a' and 'a — +' with the arrow pointing towards the respective center. 
We see that the central MT-sphere is so large, that it overlaps the screening 
sphere of the neighboring atom. Like in Fig. ^, the full curve shows the MTO 
(<^°), and the dot-dashed ((^"F), the dotted (<p°°l"), and the dashed (■0°) 
curves show the three terms in the renormalized version of equation (^5|). 
The dot-dashed and the dotted curves are identical with those in Fig. ||, ex- 
cept for the normalization; they are the outwards- integrated solution {ip°'Y) of 
the radial Schrodinger equation, continued by the inwards-integrated solution 
[(p° °-Y) for the flat potential. These two curves have been deleted outside the 
central MT-sphere where their contribution to the MTO (^) cancels. The 
inwards integration ends at the screening sphere, inside which ip° ° -with 
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j'^ truncated- cancels its own-part, ni (nr) /ni (na) , of the screened spherical 
wave, ip, shown by the dashed curve (see Eqs. ( |36| ) and ji^)). Neither of these 
cancelling parts are shown in the figure, and the dashed curve inside the cen- 
tral screening sphere therefore merely shows the contribution to the screened 
spherical wave from the inactive channels > 3). Due to the j^-truncations, 
the screened spherical wave has kinks at all screening spheres and, inside 
these spheres, only the contribution from the inactive partial waves -which 
are regular solutions of the radial Schrodinger equations- remain. The full 
curve is the MTO, which is identical with the screened spherical wave outside 
its own MT-sphere. At its own screening sphere, its kink differs from that of 
the screened spherical wave due to the truncation of the ^^-contribution to 
Compared with the bare MTO in Fig. |, the screened MTO in Fig. | is 
considerably more localized, even though a negative energy was chosen. 

If one demands that the valence band -as well as the lower part of the 
conduction band- of Si be described from first principles using merely the 
minimal 4 orbitals per atom, one cannot use a set with p orbitals such as 
those shown in Figs. ^ and |[ the d-MTOs must be folded into the envelopes 
of the remaining sp set by use of the appropriate structure matrix obtained 
from Eq. ( ^2|) with the choice ( ^3|) for the Si d-channels. The corresponding 
Si piii-MTO is shown in Fig. ^. Little is changed inside the central screening 
sphere, but the tail extending into the nearest-neighbor atom has attained 
a lot of c?-cliaracter around that site, and the MTO is correspondingly more 
delocalized. 

The Si piii-MTO for use in an sp MTO basis constructed from the con- 
ventional Si-|-E potential -for which the diamond structure is packed bcc 
with equally large space-filling spheres- is obtained by down-folding of the 
Si d and all empty-sphere channels. It turns out to be so similar to the one 
obtained from the fat Si-centered potential shown in Fig. U, that we will not 
take the space to show it. 

Si p kinked partial wave 



0.4 



0.2 



0.0 

Si <-a a^s Si <-a 
Fig. 4. Si Pill member of a screened minimal sp-set of Oth-order MTOs (see text). 
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Whereas the bare MTO in Fig. is what has always been called a bare 
MTO, the screened ones in Figs. H^and ^ look more like a partial wave, 
ipY, with a tail attached at its own screening sphere -and with kinks at all 
screening spheres. Hence the name 'kinked partial wave' given in Ref. In 
this original derivation, kinked partial waves with a = s < t were considered 
first, and only later, the limiting case a — > gave rise to a painful exercise. 
The kinked partial waves have in common with Slater's original Augmented 
Plane Waves (APWs) |59j, that they are partial waves, (e, r) Y, of the proper 
energy inside non-overlapping spheres, which are joined continuously -but 
with kinkS" to wave-equation solutions in the interstitial. In that region, the 
APW is a wave-equation solution with a given wave-vector, whereas the MTO 
is a solution with the same energy. Moreover, whereas the APW method uses 
identical potential and augmentation spheres, this is not the case for MTOs. 

If -for the third time in this section- we make a linear combination of 
MTOs -this time defined with kinks- and demand that it solves Schrodinger's 
equation, then the condition is, that the kinks -rather than the j'"-functions- 
from the tails should cancel the ones in the head. This condition is of course 
equivalent with the one for ^"-cancellation. Nevertheless, let us express the 
KKR equations in this language because it will turn out to have three further 
advantages: The artificial dependence on k = y/e and the associated change 
between Neumann and decaying Hankel functions will disappear, there will 
be a simple expression for the integral of the product of two MTOs, and 
we will be led to a contracted Green function of great use in the following 
section. 

Since the kinks arise because the j"-functions are truncated inside their 
screening spheres, the kink in a certain active channel of an MTO is propor- 
tional to the slope of the corresponding j "-function at a+. An expression for 
this slope is most easily found from the Wronskian, which in general is de- 
fined as: [/ (r) g' (r) — g (r) /' (r)] = {/, g}^ , and is independent of r when 
the two functions considered are solutions of the same linear, second-order 
differential equation. As a consequence, {n,j°'} = {n, j — ntana} = {n,j} = 
—K~^, and therefore: 

or (e, r) /dr\^^ - - [a'nn (na)] . (47) 

We now define the elements K'^,j^, j^j^ (e) -where R'L' and RL both refer 
to active channels- of a kink matrix |ll9| , ^ as a^/^/ times the kink in the 
i?'L'-channel of 0^^(e,r/j). From the expression for dj"' /dr\^^, the last 
forms of the spherical-harmonics expansions (^6|) and (|3^), the definition 
( pO[ ) of the screened KKR matrix, and the renormalization (|4^), this is seen 
to be: 
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Note that this is the kink matrix as defined in Ref. [£0[, whereas the one 
defined in Ref. jl^ has the opposite sign. As presentiy defined, the energy 
derivative of the kink matrix is positive definite, as we shall se in the next 
section. 

Screening and the definition ( p^ ) have removed the spurious energy depen- 
dencies of K"^^^ (e). To see this more clearly, let us use the first -rather than 
the last- forms of the spherical-harmonics expansions (|3^) and (^9|) , which are 
also more closely related to the definition ( ^5[ ) of the MTO, and to Figs, ^and 
^: The kink matrix for ip'\ {£,vii) is — [ktlii (ko^')] ^ B%j^ (e) [nni (/coa)] 
Moreover, (e, rjj) contains the diverging term n (nr) jn (kc) in its own 
channel, which in the MTO is being cancelled by a term from 93° " (see the 
third equation ( |36| ) and (|3^)). The kink matrix for the MTO set is now seen 
to equal the one for the set of screened spherical waves, plus -in the diagonal- 
the kink in the function Lp°- — °- ^ n (nr) jn (na) . Since ip — ip° is smooth, 
this kink is the one between the radial functions 1^9° ° (e, r) and n (nr) jn {na) . 
We thus arrive at the expression: 

K%,L',RL ie) = Br'L' RL (g) ^ (49) 

+aRL [d {ni (e, a)} - d {ip° (e, a)}] Sr>rSl>l 
gaO {if'^ (e, a)} 6wb6l'l 



d 



= Bb'l',rl{£) - aBLd{ip° {e,a)} Sb'rSl'l, (50) 

in terms of the logarithmic-derivative function at the screening sphere of 
the inwards-integrated radial function, d {(^^ (e, a)} = d\n (e, r)\ /d\nr\^. 
Remember that RL and R'L' refer to active channels. 

In the third line of (|4^) we have pointed to the fact that the first, potential- 
independent part of the kink matrix is a\, times the outwards slope of 
the screened spherical wave and in ( ]50| ) we have denoted this slope matrix 
^R'L' RL (^) ■ l^ote that, as presently defined, this slope matrix is Hermi- 
tian and equals a^'L' times the non-Hermitian slope matrix defined in Refs. 
[ p^|j20| ; moreover, the transformation from to B"^ is not quite (^), but 
differs from it by the term a91n \ni {Kr) \ /91nr|^. We may switch from Neu- 
mann to Bessel functions, using again that ji (na) = ni (Ka)tana, and that 
{ji"} — 1/^- We get: 

, , tan a (e) r „^ / N ^ n tan a(e) ^ r , 
B (e = 7^ [B" (e) - k cot a (e) — ^ + ad {j (na)} 

= . / , [B (e) -\- KCota(e)V^ . / , -f ad{j (na)} , (51) 
J (ko) j (Ka) 



where the last equation has been obtained with the help of (|32|), and where 
B (e) = B"^" (e) is the bare KKR structure matrix (|2|). The matrix B (e) -\- 
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K cot a (e) is the bare KKR matrix for the background-potential and has 
dimension (A + /) ; it only truncates when aj (e) = rji (e) — 0, as it happens 
for high I. 



Computational Procedure. The recipe for a computation could be: Solve 
the radial Schrodinger equations outwards, and then inwards to a ~ 0.8i, 
for all channels up ^ < 3. Then, compute the Green matrix of the back- 
ground, G"^" (e) = [_B"^" (e) -I- n cot a (e)] , by inversion in real space, 
choosing the strong screening just mentioned, i.e. nearly touching screening 
spheres for all spd{f) channels. This gives the strongly screened structure 
matrix, (e) or B'^ (e) , according to (|5l|), and the KKR matrix, K'^ (e) or 
-fC" (e) , for the real potential in the strongly screened representation accord- 



ing to (|40| ) or (|50| ). For a crystal, Bloch-sum the KKR matrix. Now, invert 
this matrix in real space to obtain the Green matrix, G" (e) = K" (e)" or 
G" (e) = K"^ (e) . Next, choose the physically and chemically motivated 
screening (/3) and rescreen the Green matrix to the downfolded representa- 
tion, (e) or (e) , using the scahng relations js^ ) or derived be- 
low. As will be explained in the following Sect. |[ this should be done for 
a number of energies. In addition, one will need the first energy derivatives 
& (e) . The latter may be obtained from K'^ (e) via numerical differentiation 
of the weakly energy dependent structure matrix, (e) , and calculation of 

if"' (e,r) r'^dr — ip]^^ (Sj'^) r^dr for the energy derivative of the loga- 
rithmic derivative function in (^), as will be shown in (|60|)-(|62|) below. With 
this iC^ (e) , compute G" (e) from ( |6^ and, finally, rescreen to G^ (e) using 
the energy derivative of given below. 

In order to evaluate the wave function (^), one needs in addition to 
^A'A i^) J block Bjj^ (e) , and this may be obtained from (^). 

The relation of the screening constants, the structure matrix, and the 
KKR matrix to those -see and (^0|)- of the conventional ASA is sim- 
ple, but not as straightforward as the a-to-a transformations of the present 
section, so for this topic we refer to Refs. p^ , pOt . 

This completes our exact transformation of the original KKR matrix ( p5[ ) 
which has long range and strong energy dependence (e, k) has poles at 
the free-electron parabola: X^gI^"*"^! ^ screened and renormal- 

ized KKR matrix which -depending on the screening- may be short ranged 
and weakly energy dependent. The kink matrix is expressed in terms of a 
slope matrix, which only depends on the energy and the structure of the 
background, and the logarithmic derivatives of the active radial functions 
extrapolated inwards to the appropriate screening radius. 



2.4 Re-screening the Green matrix 

In the ASA, it is simpler to re-screen the Green matrix ( p^ ) than the structure 
matrix M) , because the former involves additions to the diagonal and energy- 
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dependent rescaling of rows and columns, but no matrix inversions. The same 
holds for the fully energy-dependent matrices of the 3rd-generation, as may 
be seen from (32) or ( [5l] ) for the structure matrix. For the Green matrix (^0|), 
we get with the help of ( pT] ) and a bit of algebra: 

G" (e) EE if" (e)"^ = K-4ana(e) [1 - tan a (e) cot 77 (e)] 

+ [1 - tan a (e) cot r] (e)] G"=° (e) [1 - tan a (e) cot 77 (e)] , 

which has the form (0). Solving for G"^" (e) and setting the result equal to 
G^ (e) yields the following relation for re-screening of the Green matrix: 

tan77'^(e)^^ tan 77'' (e) tana (e) — tan/3 (e) tan?7'' (e) 
~ tan 77" (e) tan 77" (e) k tan 77" (e) ' 

(52) 

In a-language, where according to (^8|): G" (e) = — an (na) 0°" (e) nn (na) , 
the diagonal matrices in ( p2[ ) become [71 (k&) /7i (ko)] [tan 77'' (e) / tan 77" (e)] 
and Kn (kq) n (Kb) [tana (e) — tan/3 (e)] and may, in fact, be expressed more 
simply in terms of the inwards-integrated radial wave function, renormalized 
according to (^6|). In order to see this, we first use the form (^Sj): 

ip°''{s r) = "(^'^) tan?7° (e) - j" {e,r) 
' n (ku) tanry" (e) 

and then evaluate this at the screening-radius b : 

71 (k6) tan?7" (e) — j" (e, &) n (Kb) tan rj^ (e) 
' 71 (kc) tan77" (e) 71 (ko) tan 77" (e) 

To obtain this result, we have also used: 

j" (e, 6) = j (Kb) — 77 (nb) tan a (e) = 77 (k6) [tan (3 (e) — tan a (e)] , 

from ( pO| ) and (^2|). The second, readily computable function is that solution 
of the radial wave equation which vanishes at a with slope : 



Q_2Qja (-g^ j,^ /dr\^ 



r{e,r)^ ^~nn{na)r{s,r). (53) 



Evaluation at r = 6 yields: 

(e, 6) = — K77 (ko) j" (e, &) = kti (ko) 77 (k&) [tana (e) — tan/3 (e)] , 



which is the second function needed. Hence, we have found the following 
simple and practical scaling relation for re-screening of the Green matrix: 



G" (e) = ^° (£, b) is) ip° (e, b) + f (e, b) ^° (e, 6) 



(54) 
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2.5 Green functions, matrix elements, and charge density 

The kinked partial wave is the solution of the inhoniogeneous Schrodinger 
equation: 

{n - e) rwL' (£, r) = - ^ <5 {rn - a^x) >l (ffl.) i^Sx.fl'L' (e) , (55) 

RL 

provided that we define the MTO (|3^) the 3-fold way indicated in Figs. ||-^, 
and therefore -for the MT-Hamiltonian H (||)- use the radial Schrodinger 
equation ^ channel- wise. 

The kinks of the MTO are given correctly by (^5|) , but the proper MTO 
does not solve Schrodinger's differential equation in the shells between the 
screening and the MT-spheres; here we need the 3-fold way. This way must 
not be an approximation: For instance, when applied to those linear combi- 
nations of MTOs which solve the KKR equations -and hence Schrodinger's 
equation- equation (^5|) is correct (and yields zero), because for each ac- 
tive channel. A', the two solutions, Va' J2a (^j ^r) ^^"^ Va' (^^ ^R') ^a'^ 
of the radial wave equation match in value and slope at ajii , and there- 
fore cancel throughout the shell sri — an'L'- Expressed in another way: For 
energy-dependent MTOs, kink-cancellation leads to cancellation of the triple- 
valuedness. For the energy-independent NMTOs to be derived in the next 
section, special considerations will be necessary. 

Solving (H) for 5 {rR - orl) Yl {vr) , leads to: 

{n~e)Y, ^R'L' (e, r) G^^.^r^ (e) = -5 (tr - orl) Yl [vr) (56) 

R'L' 

which shows that the linear combinations 

IRL (e, r) = ^ cj>%L' (e, r) G%l',rl (e) , (57) 

R'L' 

of MTOs -all with the same energy and screening- is a contraction of r' onto 
the screening spheres (r' — > orl, RL) of the Green function defined by: 

(Hr - e) G (e; r, r') = -(5 (r - r') . 

The contracted Green function ^r^ (^i '') kink 1 in its own channel and 
kink in all other active channels (7^ RL) . This function is therefore a solu- 
tion of the Schrodinger equation (defined the 3-fold way) which is smooth ev- 
erywhere except at its own screening sphere. 7^^ (e, r) is usually delocalized, 
and when the energy, e, coincides with a pole, Sj, of the Green matrix, 
7^^ (£,r) diverges everywhere in space. This means, that when e — then 
the renormalized function is smooth also at its own sphere, and it therefore 
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solves Schrodinger's equation. In vector-matrix notation, equations (|55|) and 
( p6| ) become: 

{n~e)cj,^ {e,r)^-d^ {r) K^e) , 
{n - e) r {s, r) G"^ (e) ^{H-e) ^ (e, r) - -S^ (r) , 

where we have defined a set of spherical harmonics on the a-shells with the 
following members: 

d]fj^{rH) = d{rR~aRL)YLirR). (58) 

If expressed in real space, our Green matrix, (e) , is what in multiple- 
scattering theory is usually called the scattering path operator and de- 
noted T (e) . In the 2nd-generation LMTO formalism, it was denoted g (e) , 
but in the present paper we denote matrices by capitals. 

Since in the 3-fold way, an MTO takes the value one at its own screening 
sphere and zero at all other screening spheres, expression ( |55| ) yields for the 
matrix element of 7i — e with another, or the same, MTO in the set: 



(e) \n - e\ ^RL (£)> = -K^R'L'ML (e) = -G^^^^r^^ {e)-' , (59) 



which says that the negative of the kink matrix is the Hamiltonian matrix, 
minus the energy, in the basis of energy-dependent Oth-order MTOs. 

For the overlap integral between screened spherical waves, with possibly 
different energies and in the interstitial between the screening spheres, defined 
channel-by-channel, we obtain the simple expression [ p^ : 

(rR-.ieVrRAe)) = """"-"^^'^llf""'""^'^ (60) 

by use of Green's second theorem, together with expression (|5^) for the sur- 
face integrals. Note that, neither active channels different from the eigen- 
channels, i?'L' and RL, nor the inactive channels contribute to the surface 
integrals. The reasons are that iPr'l> (e', r) and V'^x (^i vanish on all 'other' 
screening spheres, and that they are regular in the inactive channels. The 
latter means that, in the inactive channels, the 'screening-sphere interstitial' 
extends all the way to the sites (a/ 0). For the overlap integral between 
kinked partial waves, the 3-fold way yields: 

irWL' (£') I ^RL (£)> ^ irR'L' (£') I ^RL (£)> + Sr,rSl,l X 

r (£', r) ip^L (£, r) r'dr - H {s' , r) (e, r) r'dr) 

■ } ■ — ' Kr'L'.rl (e) if e ^ £• (61) 

e — e 
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For the overlap matrix for the set of contracted Green functions, this gives: 
(7- is') I r {e}) - ^9lS5l_^lS5l (62) 
^ -G"(e) = (e) ii:" (e) G" (e) if e' ^ e. 

We see that S° (e) , (e) , and G" (s) are Hermitian, just like B'' (e) , K" (e) , 
and G" (e). Whereas (e) and iiT" (e) are positive definite matrices, that 
is, their eigenvalues are positive or zero, G"^ (e) is negative definite. For well- 
screened MTOs, the logarithmic derivative functions in the diagonal of the 
kink matrix ( |50| ) depend more strongly on energy than the slope matrix. The 
way to compute the energy derivative K"^ (e) is therefore to compute 5° (e) 
by numerical differentiation, and the remaining terms by integration as in 

(ill)- 



In the following we shall stay with the normalization (|4J)-(46) denoted by 
Latin -rather than Greek- superscripts and shall rarely change the screening. 
We therefore usually drop the superscript a altogether. Some well-screened 
representation is usually what we have in mind, but also heavily down-folded 
-and therefore long-ranged- representations will be considered. In those cases, 
some parts of the computation must of course be performed in the Bloch -or 
k-space- representation. 

The wave function is Wi (r) = (ci, r) Cj , where the eigen(column)vector 
Ci solves the KKR equations, K (si) Ci = 0, and is normalized according to: 
1 = c\k {si) Ci, in order that (S'^ | S'^) = 1. From the definition ( p5| ) of the 
MTO, we see that an accurate approximation for the charge density, which 
is consistent with the 3-fold way and, hence, with the normalization, has the 
simple form: 

p (r) = (r) + ^ [pI (r«) - p^ (r^)] (63) 

R. 

where the global contribution is: 

p'^ (r) = ^ ^ r i^RL (e, vr) Frl.r'L' (e) V'fl'L' (e, vr,)* de (64) 

RR' LL' •' 

and the local contributions, p'^ (r^) — p^ (^i?) , which vanish smoothly at 
their respective MT-sphere, are given by: 

ipRi (£, r) Frl^rl' (e) VRV (e, r) de 
PR (r) = E ) ^i' ) r ^Ri ^RL,RL' is) H>°Rv (£, r) de . (65) 

LL' •' 
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The common density-of-states matrix in these equations is: 

occ ^ 

rRL.R'L' (e) = CRL,iS (e - Ej) c^/ i = -ImGi^L^ffL' (e + iS) ■ (66) 

The approximations inherent in ( |6^ ) are that all cross-terms between prod- 
ucts of ■0-, (^-, and (^"-functions, and between (p- or (/3°-functions on different 
sites are neglected. 



3 Polynomial MTO approximations 

In this section we shall show how energy- mdependent basis sets may be 
derived from the kinked partial waves, that is, how we get rid of the energy 
dependence of the MTOs. Specifically, we shall preview the generalization 
[ pTp^ of the 3rd-generation LMTO method mentioned in connection 

with Fig. |l|. This generalization is to an 'N'MTO method in which the basis 
set consists of energy- mdependent NMTOs, 

N 

X^RL (r) = E E ^R'^' I^R'l,RL;n - (67) 

n=0 R/L' 

N 

where E LR'h.,RL;n = ^R'R^L'l, 

constructed as linear combinations of the kinked partial waves at a mesh of 
A^-|-l energies, in such a way that the NMTO basis can describe the solutions, 
]Pi (r) , of Schrodinger's equation correctly to within an error proportional to 
{si — £o) {si — £i) ■•• (£i — £n) ■ Note the difference between one-electron en- 
ergies denoted and e^, and mesh points denoted £„ and Em, with n and m 
taking integer values. The set, x^^"°'' (r) , is therefore simply 0(eo,r) , and 
this is the reason why, right at the beginning of the previous section, (e, r) 
was named the set of Oth-order energy- dependent MTOs. For iV > 0, the 
NMTOs are smooth and their triple- valuedness decreases with increasing N. 
For the mesh condensing to one energy, e,j, the NMTO basis is of course con- 
structed as linear combinations of (p (e,^, r) and its first N energy derivatives 
at £1,. For 7V=1, this is the well-known LMTO set. 

The immediate practical use of this new development is to widen and 
sharpen the energy window inside which the method gives good wave func- 
tions, without increasing the size of the basis set. One may even decrease the 
size of the basis through downfolding, and still maintain an acceptable energy 
window by increasing the order of the basis set. The prize for increasing N 
is: More computation and increased range of the basis functions. 
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3.1 Energy-independent NMTOs 

What we have done in the previous sections -one might say- is to factorize out 
of the contracted Green function, 7 (e, r) , some spatial functions, (jjRL (£, r) , 
which arc so localized that, for two energies inside the energy-window of 
interest, the corresponding functions, and 0HL(£',r), cannot be 

orthogonal. In other words: The kinked partial waves are so well separated 
through localization and angular symmetry that we need only one radial 
quantum number for each function. 

Now, we want to get rid of the kinks and to reduce the triple-valuedness 
and the energy dependence of each kinked partial wave retaining its RL- 
character- to a point where the triple-valuedness and the energy-dependence 
may both be neglected. This we do, first by passing from the set (l>{s,r) 
to a set of so-called iVth-order energy-dependent MTOs, X^''^'' {^^ ^) ) whose 
contracted Green function, 

AT 

x(^) (£, r) G (e) ^ (/. (£, r)G{e) (£„, r) G {e) , (68) 

differs from (p (e, r) G {e) by a function which remains in the Hilbert space 
spanned by the set d (e, r) with energies inside the window of interest, and 
which is analytical in energy. The two contracted Green functions thus have 
the same poles, and both energy-dependent basis sets, (e, r) and x'^-* (e. r) , 
can therefore yield the exact Schrodinger-equation solutions. The analyti- 
cal functions of energy we wish to determine in such a way that x'^^ (e,r) 
takes the same value, x^''^^ (r) , at the -|- 1 points, £0, ....En- With the set 
X^''^^ (e, r) defined that way, we can finally neglect its energy dependence, 
and the resulting x*'^'' (i") is then the set of A''th-order energy- independent 
MTOs. 

Other choices for the analytical functions of energy, involving for instance 
complex energies or Chebyshev polynomials, await their exploration. 

One solution with the property that Xb2 takes the same value for 
e at any of the A'' + 1 mesh points, is of course given by the polynomial: 

N 

„ , En Em 

of A^th degree. But this solution is useless, because it yields: x^^-* (r) — 0. If, 
instead, we try a polynomial of (N — l)st degree for the analytical function, 
then wc can write down the corresponding expression for the set x^^^ (i") 
without explicitly solving for the (iV -|- 1)'^ matrices a[^^ {sm) , and then 
prove afterwards that each basis function has its triple-valuedness reduced 
consistently with the remaining error oc (sj — £0) (^i — £1) ••• (e* — en) of the 
set. 
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Since we want x*-^^ (£„, r) to be independent of n for < n < A^, all its 
divided differences on the mesh -up to and including the divided difference 
of order N- vanish, with the exception of the 0th divided difference, which 
is x''^'' (i")- A-S a consequence, the A^th divided difference of x^^^ (£: r) G (e) 
on the left-hand side of ( |6^ ) is x'^^'' (r) times the iVth divided difference of 
the Green matrix. Now, the A^th divided difference of the last term on the 
right-hand side vanishes, because it is a polynomial of order N — 1, and as a 
consequence, 

^ ^' A[O...N] yA[O...N]J ■ ^ ' 

This basically solves the problem of finding the energy-independent NM- 
TOs! What remains, is to factorize the divided difference of the product 
(j){e,r)G{e) into spatial functions, (/>(e„,r), which are vectors in RL, and 
matrices, G(£„), with n — 0, ...,iV. Equivalently, we could use a binomial 
divided-difference series in terms of (p (£o: r) and its first N divided differences 
on the mesh together with G (en) and its corresponding divided differences. 

For a condensed energy mesh, defined by: e„ — + En for < n < A^, the 
A^th divided difference becomes times the A^th derivative: 



= f[O...N 



A[O...N] ' ' "' ' Nl dE^ 



(70) 



but since a discrete mesh with arbitrarily spaced points is much more powerful 
in the present case where the time-consuming part of the computation is the 
evaluation of the Green matrix (and its first energy derivative for use in Eq. 
(|62[)) at the energy points, we shall proceed using the language appropriate 
for a discrete mesh. In ( |70| ) we have introduced the form / [O...A^] because it 
may -more easily than A^^ f /A [O...A^]- be modified to include another kind 
of divided differences, the so-called Hermite divided differences, which we 
shall meet later. 

Readers interested in the details of the discrete formalism are referred 
to the Appendix where we review relevant parts of the classical theory of 
polynomial approximation, and derive formulae indispensable for the NMTO 
formalism for discrete meshes. Readers merely interested in an overview, 
may be satisfied with the formalism as applied to a condensed mesh and 
for this, they merely need the translation (|70| ) together with the divided- 
difference form of the NMTO to be described in the following. Details about 
the Lagrange form may be ignored. 

Lagrange form. We first use the Lagrange form (|148| ) of the divided differ- 
ence to factorize the energy-independent NMTO (|6S|) and obtain: 

xW (r) = G [0..N]-' , (71) 
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Fig. 5. Si P ill member of the sp-set of 0th (dottet) and Ith-order MTOs (see text 
and Eq.@). 



Here and in the following, 0„ (r) = cj) (£„, r) and G„ = G (£„) ■ Eq. ( [h] ) has 
the form ( |67| ) and we see, that the weight with which the MTO set at e„ 
enters the NMTO set, is: 

4"^^ = -wv ^" _ G [O..Nr' . (72) 



By application of ( |148| ) to the Green matrix, we may verify that these La- 
grange weights sum up to the unit matrix. For this reason, the RL characters 
of the NMTO basis functions will correspond to those of the kinked partial 
waves. 

As an example, for iV=l we get the so-called chord-LMTO: 

X^i' (r) = ^0 (r) Go (Go - G.y' + 0i (r) Gi (Gi - Go)"' 

= 00 (r) (i^i - Kay' Ki + cj)^ (r) (Ko - K^y' Kq (73) 
= 0o(r)-0([Ol],r)X[Ol]-'Xo 
^ 0(r) -^{r)k-^K. 

In this case, there is only one energy difference, Eq — £i, so it cancels out. 
In the 3rd line, we have reordered the terms in such a way that the Newton 
form, to be derived for general N in (|8^) and ( p9| ) below, is obtained. In 
the 4th line, we have condensed the mesh onto whereby the well-known 
tangent-LMTO is obtained. The latter is shown by the full curve in 

Fig.|| for the case of the Si pm-orbital belonging to an sp set. The dashed 
curve is the corresponding kinked partial wave, (j) (r) , shown by the full curve 



38 O.K. Andersen et al. 



in Fig. ^ Compared to the latter, x'-^'' (r) is smooth, but has longer range. 
The strong contributions to the tail of the LMTO from (r) 's on the nearest 
neighbor are evident. It is also clear, that for computations involving wave 
functions -e.g. of the charge density- the building blocks will rarely be the 
NMTOs, but the kinked partial waves, 0„ (r) , which are more compact. 

One might fear that the discrete NMTO scheme would fail when one of 
the mesh points is close to a one-electron energy, that is, to a pole of the 
Green matrix, but that does not happen: If one of the G„'s diverges, this 
just means that the corresponding Lagrange weight is 1, and the others 0. 
Hence, in this case the NMTO is just (/>„ (r) , and this is the correct result. 
Moreover, the kink of this single 0„ (r) does not matter, because in this case 
where G (e) is at a pole, the determinant of its inverse vanishes, so that the 
kink-cancellation equations, KnCn — 0, have a non-zero solution, c„, which 
yields a smooth linear combination, </>„ (r) c„, of NMTOs. 



Kinks and triple- valuedness. The energy- independent NMTOs have been 
defined through ( |68| ) and ( |69| ) in such a way that x*"^' (^i i") ~ X^^-' (r) oc 
(e — So) ... (e — Sn)- We now show, that also the kink-and-triple-valuedness 
of x*'^-' (r) is of that order, and therefore negligible. 

The result of projecting the energy-dependent MTO onto Y^/ (f_R') for an 
active channel was given in (|3^) for its own channel, and in (^) for any other 
active channel. Together, these results may be expressed as: 

[k cot ?7gi (e) Sr'rSl'l + B%^, ,^^ (e)] 
or, in terms of the renormalizcd functions (|4^), (|45|), (|4^), and (^), as well 



as the kink matrix defined in (48), as 



Here, like in (|3^ ) and (|39|), contributions from MT-overlaps -which are ir- 
relevant for the present discussion- have been neglected. Without kinks and 
triple- valuedness, Vr/lxPrl (^i ^r) would be given by the first term, and the 
kinks and the triple-valuedness are therefore given by the second term: 

Tn,L'rRL (e, rn) - j^^l' (e, ^fl/) ^S'L',m (e) • (74) 

This vanishes for those linear combinations of MTOs which solve the kink- 
cancellation conditions. 

What now happens for the energy-independent approximation, x*-"-* (r) = 
(/>o (r) , to the Oth-order energy-dependent MTO, x^°^ (Ej r) = (/) {e, r) , is that 
the former has kinks and triple-valuedness, but both are proportional to 
K (sq) which -according to (|5^)- is proportional to Ti — Eq and, hence, to 
Ei — Sq. The kinks and triple-valuedness are thus of the same order as the 
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error of x*^*^' (r) . Similarly, for > 0, the fact that the An' (£)'s are polyno- 
mials of {N — l)st degree, reduces the triple-valuedness of x^^^ (r) to being 
proportional to {e — Eq) ... {e — En) , as we shall now see: Multiplication of 
(H) with G" (e) from the right yields: Tcj)'' (e, r) (e) = (e, r) , and for 
the kinks and the triple-valuedness of the contracted Green function (|68| ) we 
therefore get: 

N 

Tx^^) (e, r) G (e) = f (e, r) - ^ (£„, r) A^^^ (e) . 

n=0 

Taking again the A^th divided difference for the mesh on which x'^^* (S;!") is 
constant yields: 

rx(^) (r) = r ([0...iV] ,r) [0...N]-' (75) 
= -J- ([0...iV] , r) - eo) (i?^^) - £i) ... (i?^^^ - e^v) , 

for the kinks and the triple-valuedness of the energy- independent NMTO. In 
the last line, we have used an expression -which will be proved in (|8^)- for 
the inverse of the A^th divided difference of the Green matrix in terms of the 
product of energy matrices to be defined in (pO[). At present, it suffices to 
note that differentiation of the Green function, 

G(s)^E7^' (76) 
for a model with one, normalized orbital yields: 

"j_ d^G (e) 

N\ dE^ ^_ 



E 



\N+1 



(£i - £i/) 



TV -1-1 



where the last approximation holds when the mesh is closer to the one- 
electron energy of interest, Et, than to any other one-electron energy, £j ^ Ei. 
Note that j -and not n- denotes the radial quantum number. Similarly, this 
model Green function has a divided difference on a discrete mesh of A^+I 
points, whose inverse is: 

G[O..Ar; 



E 



AT 



(77) 



ji=0 



as proved in Eq. (158) of the Appendix. We have thus seen that the triple- 
valuedness is of the same order as the error present in x^^^ (r) due to the 
neglect of the energy-dependence of x''^'' (Sji") . 
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The radial function (e, r) in ( fz^ ) vanishes for r < a, where it has a kink 
of value and it solves the radial wave equation for r > a. As shown in 

, its expansion in powers of r — a > is: 



rf (e, r) 



a 

1 
5! 



/(/ + !) (r 



3! 



18; (; + l) + {l + l)~ea'^y 
This means the A^th divided-difference function entering (|75|) satisfies: 

The fcmfc and triple-valuedness (|75| ) in the s — a shell of x^^'' (r) is thus 
proportional to (r — a)^^^^ H^lLo ~ ' ^'^'^ -'^'-'^ ^^^^ reason the energy- 
window widens as s — a decreases, that is, as the screening increases. 



Transfer matrices and correspondence with Lagrange interpolation. 

We need to work out the effect of the Haniiltonian on the NMTO set. Since the 
NMTOs with A^ > are smooth, the contributions from the delta-function 
on the right-hand side of (pq ) for the contracted Green function will cancel 
in the end. Operation on (pSj) therefore yields: 



n 



< (£, r) - x'''> is, r) G (e) = </> (e, r) sG (e) - T^x^^^ (e, r) G (e) 



J2 ^0n(r)£„G„AW(e) 

^ — ^n— 



and by taking the A^th divided difference for the mesh on which x'"^'' (Sii") 
is constant, we obtain: 



H7([0...A^],r) = T^x'^^ (r)G[O...A^] = {(t>eG) {[0...N] ,r) 
= 7([0..A^-l],r) + eNl{[O...N],r), 



(78) 



using (150) with the choice of the last point on the mesh. Solving for the 
NMTOs yields: 



1) (r) 



(79) 



where x^^ (r) = 7 ([O..A^ — 1] , r) G [0..N — 1] ^ is the energy-independent 
MTO of order A" — 1, obtained by not using the last point. Moreover, 



E^^^ = EN + G[O..A^- l]G[O...Ar]~^ = (eG) [O...A^] G[O...A^]"^ 

N ^ N 



SnGfi 



G[O...A^]" 



(80) 



n=0 
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is the energy matrix which -in contrast to (r)- is independent of which 

point on the mesh is omitted. The first equation (^) shows how to compute 
E^-'^'^ and the last equation shows that iJ^^^ is the energy weighted on the 
O...A^-mesh by the Lagrange matrices ([f2[). For a condensed mesh, the results 
is the simple one ( p^ quoted in the Overview. 

We now consider a sequence of energy meshes, starting with the single- 
point mesh, Sq, then adding £i in order to obtain the two-point mesh Eq, Si, 
then adding £2 obtaining the three-point mesh eg, £1, £2, a.s.o. Associated 
with these meshes we obtain a sequence of NMTO sets: the kinked-partial 
wave set, x^°^(r)>the LMTO set, x^^UO^the QMTO set, x^^Ui") , a.s.o. 
Working downwards^ we thus always delete the point with the highest index. 
Equation (|7^) now shows that Ti. — em may be viewed as the step-down 
operator and E''^^^ ~ en the corresponding transfer matrix with respect 
to the order of the NMTO set. 

In this sequence we may include the case A^=0, provided that we define: 

i;(°)-£o = -i^(£o) and x'"^^ (r) = -5(0. (81) 
-I- 1 successive step-down operations on the NMTO set thus yield: 

{n - £0) ...{n- en) X ^""^ (r) - S (r) - £0) ... - en) 

which, first of all, tells us that one has to operate N times with -that 
is, with V^^- before getting to the non-smoothness of an NMTO. This is 
consistent with the conclusion about kinks and triple-valuedness reached in 
the preceding sub-section. Secondly, it tells us that the higher the N, the more 
spread out the NMTOs; if we let r (M) denote the range of the E<^^^ -matrix, 
then the range of the NMTO is roughly J2m=o (^) • 

The product of E''^"^ — £0 and all the transfer matrices on the right- 
hand side of the above equation is seen from (^) and ( ^l| ) to be simply: 
— G [0...iV] ^ . Hence, we have found the matrix equivalent of the elementary 
relation (|77|): 

-G [O...N]-' - - £0) - £1) ... (e(^) - £^) . (82) 

The other way around: Recursive use of ( |8^ ) with increasing N, will generate 
the transfer matrices and will lead to the first equation (^). Note that al- 
though the order of the arguments in the divided difference on the left-hand 
side is irrelevant, the order of the factors on the right-hand side is not, since 
the transfer matrices do not commute. That G [0...N] is Hermitian, is not so 
obvious from ( ^4) either. Fi nally , we may note that G [0..n — l,n+ 1..N] is 
not defined by^|) but by (pi?) ): 



G[0..n-l,n+l..N] = G[0...7V-1] -I- {en - En) G [0....N] 
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Relation (|8^) now gives the following form for the Lagrange weights 

^ (£„ - £o) •• (£n - £n-l) (e™ - En+l) ■• (£n - Eat) ' 

(83) 



and this is seen to pass over to the classical expression (145) for the Lagrange 
coefficients if we substitute all energy matrices by the energy: e. This 

correspondence between -on the one side- the set (j) (e, r) and the Lagrange 



polynomial approximation (145) to its energy dependence (Fig. |^) and -on 



the other side- the set x'"^'' (r) expressed by ( |67| ) with the matrix form ( p3[ ) 
is conceptually very pleasing. What is not so obvious -but comforting- is 
that the Hilbert space spanned by the NMTO set is invariant under energy- 
dependent linear transformations, (j) (e, r) = (e, r) T (e) , of the kinked par- 
tial waves. This will be shown in a later section. 

By taking matrix elements of (fzol), the transfer matrix may be expressed 

as: 

This holds also for A^=0, provided that we take the value of x^^"^ (r) at its 
screening sphere to be (f° ° = 1 -as dictated by the 3-fold way- so that 
{x^^^ I X*'~^^) = 1- The form ( p^ shows that the transfer matrices with N > 1 
are not Hermitian, but short ranged, as one may realize by recursion starting 
from A^=0. Finally, it should be remembered that the NMTOs considered 
sofar have particular normalizations, which are not: (x^^^ I X^^^) — 1j ^'^d 
so do the transfer matrices. We shall return to this point. 



Newton form. Instead of using the Lagrange form (148) to fa ctorize the 
NMTO (|69|), we may use the divided-difference expression (149). With the 
substitutions: / (e) ~> G (e) and g (e) —> <j) (e, r) , we obtain the Newton form 
for the NMTO which most clearly exhibits the step-down property 



X^^^ (r) = <t>{[M..N] , r) G [0..M] G [Q...N]-' 

= 4>N (r) +4>{[N- 1, N] , r) - en) + .. (85) 

.. + ([0...7V] , r) - ei) .. - en) , 

since, from ( |55| ) and ([78|), 

{U - En) (I>n (r) = ~Sn,oS (r) Kq, 

{n - en) <j){[M...N] , r) = cj){[M..N - 1] , r) . (86) 

We thus realize that the energy matrices in the Newton series for the NMTO 
set are the matrices for stepping down to the sets of lower order. For some 
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purposes, the 'reversed' series, obtained from (149) with f [e) ^ (j) (e^ r) G (e) 
and g{e) ^ G (e): 



AN) 



= 0o(r)+0([Ol],r)(i?W-eo) + .. (87) 
.. + ([0...7V] , r) - .. (ii;^ - £o) , 

is more convenient. This expression clearly exhibits the correspondence with 
the Newton polynomial approximation (146) to the energy dependence of 
(j) (e, r) . Conceptually, a divided-difference series is more desirable than the 
Lagrange series, because the Lagrange weights ( ^3| ) 'fluctuate wildly' as a 
function of n, taken in the order of monotonically increasing energies. 

For a condensed mesh, ( |85| ) and (|8^) obviously reduce to one-and-the- 
same matrix-equivalent of the Taylor series for <f) (e, r) : 



xW(r) - 0(r) + 0(r)(£;W-e,) + 



1 (AT) , 

.. + —</> (r)(£;W-e, 



and 



becomes: 



(N-M) 

0(r) 



(N-M-l) 

0(r) 



(7V-M)! {N-M -1)1' 



Readers used to the LMTO-ASA method, where -according to (11)- the 
KKR matrix is basically the two-center TB Hamiltonian, may not like the 
thought of having to differentiate its inverse, the Green matrix, with respect 
to energy. (The computer seems to work well with the formalism based on 
the Green matrix). Such readers might therefore prefer an NMTO formalism 
in terms of kink matrices. For a discrete mesh many ugly relations exist, but 
the one relation which is conceptually pleasing is the following: 



= 

+ K [01] -ea)+.. + K [Q..N] (e'^^^ - e^^i) .. (i;^ - eo) , 

because it looks like the matrix form of the secular KKR equation: (e) | =0. 
This relation may be obtained by taking the A^th divided difference of the 
equation: K (e) G (e) = 1, using the binomial expression ( |149| ) for a product 
like in (|87|), but with K (e) substituted for (e, r) , and multiplying the result 
from the right by G[O...N]~^ . To find the transfer matrices from (|8^), we 
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may solve for — eq and do recursion starting from A^=l. The results are: 
E*^^^ -Eo ^ -K[Oir^ Kq -> ~k-^K, 

E^^^ -Eo = - {^K [01] + K [012] - ei) ) ^ (89) 

^ - (^k- kk-^K/2^ ^ K, 

a.s.o. These low-iV expressions are reasonably simple. For N=l, the discrete 
form is seen to be identical with ( |7^ ) and, for a condensed mesh, it reduces to 
the well-known expression for the 3rd-generation LMTO. We conclude that 
the energy matrices, E'^^^^ , are well-behaved functions of the kink matrix and 
its divided differences, up to and including Mth order. With M increasing, 
the corresponding expressions for E''^^) however become more and more com- 
plicated. The simplest expression for is therefore (|80|), the one which 
uses G-language. 



3.2 Variational NMTO method 



The NMTO set has been defined through ( |6§| ) and (69) in such a way that 
its leading errors are proportional to {e — Eq) ■■ {s — £n)- By virtue of the 
variational principle, solution of the generalized eigenvalue problem (^) with 
this basis set will therefore provide one-electron energies, Ei, with a leading 
error cx (e.; — Eq)^ .. (e^ — Eat)^ . The error of the wave function will of course 
still be of order (e.^ — Eq) .. (e^ — En) , but that is usually all right because, 
as mentioned at the beginning of the present section, the MTO scheme is 
based on the factorization: 7 (e, r) = (e, r) G (e) , where (j) (s, r) has a smooth 
energy dependence and G (e) provides the poles at the one-electron energies. 



Hamiltonian and overlap matrices. For a variational calculation, we 
need expressions for the NMTO overlap and Hamiltonian matrices, (x^^-* I X^^"*) 
and (x'^-* I'^l X^'^- From (|68|), the A^th divided difference of the contracted 
Green function ( p7| ) is: 

7^^) {[0..N] , r) = x^^^ (r) G [0..N] ^ ^" ^"^ ^ (90) 



and using now (62), we obtain for the integral over the product of the Mth 
and A^th divided differences of contracted Green functions: 

(7[0...M] I 7[0....A^]) = E E — 

n=0„'=0 Yl (£„-£,„) n (en'-^m') 
m—O.^n m' —Q^^n' 

{M+N+1) 

= -G[[O...M]..N] ^ — — . (91) 
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This is simply the negative of the {M + + l)st Hermite divided difference 
(151) of the Green matrix, as proved in Eq. (159) in the Appendix! 
Note that the meaning of a matrix equation like (p3) is: 



(7-RL (£n) I IR'L' {£n')) = -GrlM'L' [n,n'] 

= -Grl.r'l- = {iRL (£«') I 7-R'L' (e™)) ■ 

In matrix notation, that is: (7„ | 7„') — {jn' \ In) , and not: (7^ | 7„/) = 
{"fn' I 7n) • Even without the symmetry of the matrix G [n, n'] with respect 
to the exchange of n and n' , it is of course always true that 

(7-RX (£n) I IR'L' (£«')> = ilR'L' (£«') I IRL (Sn))* , 

i.e. that a matrix like (7,1 | 7„/) is Hermitian: (7,1 | 7„/) = {jn' \ 7™)^ • The 
point is, that n is an argument - not an index - of a matrix. Similarly, N and 
M are not matrix indices in (|9^) . Since the first expression ( pl[ ) is symmetric 
under exchange of N and M, because G [n, n'] is symmetric, we may choose 
M < N, and this has in fact been done in the second expression. 

From (|7|) and (|l]), we now see that the Hamiltonian matrix between the 
A^th divided differences of contracted Green functions becomes: 

(7 [0...N] \n - £jv| 7 [0...7V]) = (7 [0...iV] I 7 [0..iV - 1]) 

(2N) 

= -G[[0..iV-l]iV] ^ (92) 

Hence, we have arrived at the important results: The NMTO overlap matrix 
may be expressed in terms of the A^th-order divided difference and the (2 + 
l)st Hermite divided difference of the Green matrix as: 

(x^^Mx^^^) = -G[O...A^]"'G[[O...Af]] G[O...A^]"\ (93) 

where the -even simpler- result for a condensed mesh was quoted in the 
Overview ([l|). The Hermite derivative G [[0, A^]] is thus negative definite. 
The NMTO Hamiltonian matrix may be expressed analogously, in terms of 
a 2A^th-order Hermite divided difference: 

(Jx^^^ |W - £n\ X^^^) = -G [O.-.N]-^ G [[Q..N - 1] A^] G [O...A^]~' . (94) 

Here again, the result given in (|l^) for a condensed mesh is even simpler. 
The NMTO Green function is 

G [O...A^] {G [[Q..N -l]N]-{z- en) G [[Q...N]]]~^ G [Q.-.N] 



Expressions ( |93[ ) and ( |94| ) for the NMTO overlap and Hamiltonian ma- 
trices are not only simple and beautiful, but they also offer sweet coding and 
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speedy computation. For a crystal, and transforming to k-representation, 
one may even use the representation of contracted Green functions where the 
overlap and Hamiltonian matrices -according to (^) and (|9^)- are merely 
-G [[0...N]] and -G [[0..iV - 1] A^] . In Section | we shaU see that an energy- 
dependent linear transformation of the kinked partial waves does not change 
the Hilbert space spanned by an energy-independent NMTO set -but only 
the individual basis functions. Therefore, we might also use kinked partial 
waves (j)" (£,r) and the Green matrix G" (e) with phase-shift normalization. 

In summary: The variational NMTO scheme requires computation of the 
kink matrix and its first energy derivative at the mesh points. It delivers 
energies and wave functions which are correct to order 2N+1 and N, respec- 
tively. This lower accuracy of the wave functions is appropriate because the 
kinked partial waves are rather smooth functions of energy. For the compu- 
tation of the 9„'s entering Kn = o- (^Bn — , radial normalization- integrals 
should be used. 

As an example, for the LMTO method, the Hamiltonian and overlap 
matrices are respectively: 

xW|W-£i|x^'^) = -G[01]-^ G[[0]1] G[01]-^ 
= (£o-£i)(Go-Gi)-' (-Go + G[01]) (Go-Gi)-' (95) 

^ -G-i^G-i = -K + KK-^^k-^K, 

and 

(x^ IX^^)) = -G[01]-^ G[[01]] G[01]-^ 
= (Go-Gi)"'(-Go + 2G[01]-Gi)(Go-Gi)~' (96) 

^ = K - Kk-^^ -^k-^K + Kk-^^k-^K. 

3! 2! 2! 3! 

The result for a condensed mesh in terms of the kink matrix and its first 



three energy derivatives is seen to be almost identical to the one (|15), which 
in previous LMTO generations required the ASA. To get exactly to (1^), one 
needs to transform to the LMTO set: x^^^ (r) = x''^^ (r) -fiT"^/^, which in fact 
corresponds to a Lowdin orthonormalization of the Oth-order set. We shall 
return to this matter in Sect. ^. From the above relations we realize that 
-even for a condensed mesh and N as low as 1- G-language is far simpler 
than if- language. 



Orthonormal NMTOs. In many cases one would like to work with a repre- 
sentation of orthonormal NMTOs, which preserves the i?L-character of each 
NMTO. In order to arrive at this, we should - in the language of Lowdin - 
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perform a symmetrical orthonormalization of the NMTO set. According to 
(P3|) such a representation is obtained by the foUowing transformation: 



^(N) ^ ^(N) (j.) g [0^^^^] V_G[[0...iVr , (97) 
because it yields: 

(x^^) I x^^)) = -v/-G[[0...7V]]"'^G[[0...iV]] V-G[[0...iVr' = 1. 

Note that this means: ~G [[0..N]] = ^-G [[0..7V]] [[0-N]]. In this or- 

thonormal representation, the Hamiltonian matrix becomes 

(x(^)|H-£^|x^^)) ^ -V-G[[O...N]f'^ X (98) 

G [[0..iV - 1] TV] ^y-G[[Q...N]f\ 

To find an efficient way to compute the square root of the Hermitian, positive 
definite matrix — G [[0...iV]] may be a problem. Of course one may diagonalize 
the matrix, take the square root of the eigenvalues, and then back-transform, 
but this is time consuming. Cholesky decomposition is a better alternative, 
but that usually amounts to staying in the original representation. Lowdin 
orthogonalization works if the set is nearly orthogonal, because then the 
overlap matrix is nearly diagonal, and Lowdin's solution was to normalize 
the matrix such that it becomes 1 along the diagonal and then expand in the 
off-diagonal part, O : 

VTTo"' = i-^o + lo^ - ... (99) 

2 6 

This should work for the NMTO overlap matrix (||) when the NMTOs are 
nearly orthogonal, but it hardly works for — G [[O...A'^]] . There is therefore 
no advantage in pulling out the factor G [O...Af] , on the contrary. The other 
way around: In order to take the square root of — G [[0...iV]] , we should find 
a transformation, T, such that T't^G [[O...Af]] T is nearly diagonal, and then 
perform the Lowdin orthonormalization on the latter matrix. We shall return 
to this problem in Sect. 



One-orbital model: switching behavior of H^-'^\ and the vari- 

ational energy. Our development of the NMTO formalism has been focused 
on its matrix aspects and, through the introduction of energy matrices and 
by pointing to the correspondence with classical Lagrange and Newton in- 
terpolation of the energy-dependent kinked partial waves, we have tried to 
make the reader accept the seemingly uncomfortable fact, that the quantities 
of interest do arise by energy differentiations of a Green matrix. 
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Fig. 6. Switching beliavior of i?'^' (e^) = H^-^^ (e^) for the orthonormal one- 
orbital model defined by Eq. with 4 radial levels: Ej = 0, 1, 2, 3. 
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Let us now illustrate the Green-function aspects by considering the 1 x 
1 Green matrix (|7^) for one, normalized orbital: x''^'' (^) = (i') with 
^(JV)|2^ _ Note that in this model, j runs over the one-electron en- 
ergies, which is a different set -with much larger spacing- than the energy 
mesh whose points are denoted n and m. For a crystal, and using Bloch- 
symmetrized NMTOs and Green matrices, x''^'' (k, r) and G (e, k) , this would 
be an s-band model with j being the radial quantum number. We want the 
NMTO to describe the i-band and therefore choose the mesh between Si-i (k) 
and £,;+! (k) . In the following we shall drop the Bloch vector and not neces- 
sarily consider a crystal. 

We first demonstrate how = H^^^ -in this case a 1 x 1 Hamiltonian 
(see Sect.||)- expressed in terms of ratios of energy derivatives of a Green 
function, with its singular behavior, produces correct results for the one- 
electron energy and how, when the mesh is swept over a large energy interval, 
//(^) switches between bands with different radial quantum numbers. From 
(|0|) and (0) we get: 



En 



G[0..iV- 1] 
G[O...N] 



1 + En 



(e.j - En) 



m— Sj - 



tN 



' nm=0 Sj- 



Hence, for the model and an energy mesh with iV + 1 points, equals Si 

to order TV, with an error proportional to (e^ — Eq) .. (e^ — Eat) , which for 
a condensed mesh becomes (e^ — e^)^^^ . In Fig. ^ we show (e^) for 

= 1 to 6, computed from the above expression for a four- level model with 
Ej = 0, 1, 2, and 3, and a condensed mesh. We see that H^'^'^ (e^) behaves 
as it should: It switches from one level to the next, with the plateau around 
each level flattening out as N increases. For N odd, the switching-curve is 
step-like and, for N even, the switching is via —oo — > -\-oo. This comes from 
the ability of the denominator in the expression for H^^^ to be zero when 

-I- 1 is odd. An energy- independent orbital, as considered in the present 
model, can of course only describe one band. With the NMTO defined for 
a mesh condensed onto a chosen energy e^, we want to describe the band 
near as well as possible -also if the distance to the next band is small- 
and with a result which over a large region is insensitive to the choice of e^,. 
In a multi-orbital calculation, we should fold down those channels which are 
switching in the energy range of interest into the screened spherical waves. 
This will remove schizophrenic members of the NMTO set and prevent the 
possible occurrence of ghost bands. 

In the one-orbital model, the estimate of a true, normalized wave function, 
(ci, r) , is the A^th-order muffin-tin orbital: x*'^'' (r) = 'YlH 4'n (r) Ln^^ ■ If we 
now use ([76[) and (77) to evaluate expression (|7^) for the Lagrange weights, 
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where ^i^'' (e) is the Lagrange polynomial (145) of degree N. We have there- 
fore reached the conclusion that -in our orthonormal model, and to leading 
order- the wave function is the energy-dependent MTO, (^(e,r) , Lagrange 
interpolated over the (A^-l-l)-point mesh. 

Since the error of an NMTO set is of order A^+l, use of the variational 
principle will reduce the error of the one-electron energies, Si, from that of 
the highest transfer matrix, iJ^^^ — en, to order 2(iV-|-l). The variational 
energies are thus correct to order 2N+1. For a condensed mesh, this also 
follows trivially from (p3|)-(|9^), which show that the variational energy, with 
respect to e„, is: 



(2N) , (27V+1) 

(x^^M^-e.lx^^O _ G G 



(x^lxW) (2iV)!/ (27V + 1)! 



= H 



(2N+1) 



The odd-ordered switching curves {e^) , i/^'^^ {e^) , and H^^'^ {e^) shown 
in the left-hand panel of Fig. ^ are thus the variational estimates resulting 
from the use of respectively the 0th, 1st, and 2nd-order NMTO, that is, the 
MTO, the LMTO, and the QMTO. These curves are weU behaved. 

The expression for the variational energy in the one-band model can be 
evaluated exactly, also for a discrete mesh, and yields a transparent result. 



We use the double- mesh procedure explained in the Appendix after (151), 
and let the differences e„ = e„+Ar_|_i — e„ shrink to zero. From ( [77| ) we then 
get: 

G [[Q...N]] = -Y: /_ (100) 

1 



G[[0..7V-l]iV]--^ 



and for the variational energy (Bq): 



1 + E nLo (l: 



which of course agrees with the variational principle. 
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Fig. 7. Minimal-basis LMTO energy bands (dashed) of GaAs for two different 
choices of the screening-radii compared to the exact KKR band structure (solid). 
In the left-hand panel all screening radii were ~ 0.8t, while in the right-hand panel 
the Ga d radius was reduced to the radius of the Ga 3d core l|24|- See text. 



Treating semi-core and excited states: GaAs. An accurate description 
of the cohesive properties of GaAs requires a good band-structure calculation 
of the five Ga id}^ semi-core, the As 4s^-band, and the three As Ga 4sp'^ 
valence bands. If also the four lowest conduction bands must be described, 
one is faced with the problem of computing a band structure containing 
extremely narrow as well as wide bands over a 20 eV-region. To do this ah 
initio with a minimal Ga spd As sp basis set (13 orbitals per GaAs), has 
hitherto not been possible. 

With 1st and 2nd-generation LMTO-ASA methods one would normally 
use -RZ-dependent e^s and employ a 36-orbital-per-GaAs basis, consisting 
of the spd LMTOs centered on the Ga, the As, and the interstitial sites 
in the zincblende structure. The conduction-band errors arising from the 
choice K^=0 are so large that the combined correction is needed. Downfolding 
works for the p and d orbitals on the two interstitial spheres, but not for the 
interstitial s and the As d orbitals. With the 3rd-generation LMTO method, 
downfolding works much better, but the energy window is now screening 
dependent, and the use of i?Z-dependent £,y's is avoided because it messes up 
the formalism. 

In Fig. we show -in full lines- the exact (up to 7eV) LDA band structure 
calculated by the screened KKR method, i.e. by the 3rd-generation LMTO 
method using many energy panels and the Ga spd As sp basis. The five Ga 
Sd^" semi-core bands are at -15 eV, the As 4s^-band is around -12 eV, and the 
three As 4p^ Ga 4sp^ valence bands extend from -7 to eV. Above the gap, 
there are the four As Ap"^ Ga Asp'^ conduction bands. The dotted lines give 
results of 3rd-generation LMTO variational calculations with a condensed 
mesh and an in the middle of the three valence bands. In the left-hand 
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LMTO QMTO 




-20.0 -15.0 -10.0 -5.0 0.0 -20.0 -15.0 -10.0 -5.0 0.0 

E,(eV) E^(eV) 

Fig. 8. Mean error in eacli of the tiiree types of occupied valence bands in GaAs 
calculated with the LMTO and QMTO methods as a function of the expansion 
energy for a condensed mesh p4|. See Fig. m and text. 



figure, the screening-sphere radii for the active Ga spd and As sp channels 
were chosen at the Ga and As default values, respectively 0.82t and 0.78i, 
where t is half the nearest-neighbor distance. We see that the entire valence- 
band structure is distorted by hybridization with Ga d ghost bands. The 
dotted bands in the right-hand figure result after changing the Ga d screening- 
sphere radius to 0.35i, which is close to the actual radius of the Ga 3d core. 
Now, the band structure looks reasonable: The valence bands near Ei, are 
perfect, but the Ga 3d bands are nearly 0.5 eV to high 

That the variational LMTO method with a minimal basis and a single e^, 
cannot describe all occupied states of GaAs with sufficient accuracy, becomes 
even more obvious from the left-hand side of Fig. ^, where we show -as 
functions of e,^- the average errors of the five Ga 3d bands, those of the As 
4s band, and those of the three valence bands. The error c>c (e^ (k) — e^)^ 
of the variational energy is clearly visible for the narrow Ga 3d and As 4s 
bands. With e^^s in a narrow range around -11 eV, the variational error in 
the sum of the one-electron energies gets down to about 250 meV per GaAs. 
On the right-hand side, we show the same quantities, but obtained with the 
QMTO method. Now the errors (x (e^ (k) — Si,)^ are acceptable, and there 
is a comfortable range of Ei/'s around -10 eV where the error in the sum of 
the one-electron energies does not exceed 25 meV per GaAs. The screening- 
sphere radii chosen in these calculations were: 0.93i, 1.05i, and 0.35i for 
respectively Ga s, p, and d, and 0.89i and l.OOt for respectively As s and p. 
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Fig. 9. Like Fig. ^, but calculated using discrete meshes and as functions of the 
position of the last energy point. The first energy points were fixed at the positions 
indicated on the abscissa IMl. See text. 



In Fig. ^ we show the same kind of results, but this time obtained with the 
discrete (Lagrange) LMTO and QMTO methods. The size of the basis set, 
the screening-sphere radii, etc., were as in Fig. Hj For the LMTO method, eq 
was fixed at the position of the Ga 3d bands and the figure shows the result 
of varying the position £i of the other mesh point. The quadratic dependence 
on El of the variational energy-error oc {si (k) — Eq)^ (e^ (k) — ei)^ is clearly 
recognized. Compared with the results of the tangent LMTO method shown 
in the previous figure, those of the chord-LMTO are far superior: With ei's 
around ~5 eV, the variational error in the sum of the one-electron energies 
gets down to about 30 meV per GaAs, and yet, for N given, the method em- 
ploying a discrete mesh is computationally simpler than the one employing 
a condensed mesh. On the right-hand side of the figure, we show the QMTO 
results as functions of £2, with £0 fixed at the Ga 3d position, and £1 at the 
As 4s position. Here again, the quadratic dependence on £2 of the variational 
energy-error cx {si (k) — £0)^ (£i (k) — £1)^ {si (k) — £2)^ may be seen. We re- 
alize, that with this discrete QMTO method, meV-accuracy for the sum of 
the one-electron energies can be reached. 

Finally, in Fig |l^ we show the GaAs band structure in a wide (40 eV) 
range around the gap. Further conduction bands now appear above 7 eV and 
we needed to employ a basis consisting of the Ga spd As spdf 2E s QMTOs. 
£0 was chosen at the Ga 3d position, £1 near the gap, and £2 10 eV above 
the gap. The results of this discrete QMTO calculation shown by the dotted 
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curves agree superbly with those of a muhi-panel LMTO (=KKR) calculation 
shown in full line [^. This proves the power of the 3rd-generation NMTO 
method. 



U 

a 




Fig. 10. Energy bands of GaAs calculated with the QMTO method and the energy 
mesh indicated on the right-hand side (dashed) as compared with the exact KKR 
resuh (solid) See text. 



Massive downfolding: CaCu02. An increasingly important field of re- 
search is the electronic structure of real materials with strongly correlated 
conduction electrons. Within a given class of materials, fine-tuning of the 
interesting properties will require detailed knowledge of the single-electron 
part -the orbitals, hopping integrals and basic on-site terms- of the corre- 
lated Hamiltonian. In the previous review j2^] of the 3rd-generation 0th- and 
Ist-order differential MTO method, we demonstrated for the idealized high- 
temperature superconductor, CaCu02 with dimpled Cu02 planes, how one 
could extract low-energy, few-band Hamiltonians by massive downfolding; in 
the extreme limit: Downfolding to one Cu d^^-y^ orbital per Cu site | |2^ , p3| . 
Let us now reconsider this example in the light of the new NMTO methods. 

In Fig. |ll| the full lines in all four parts show the (same) full LDA band 
structure in a ±3 eV region around the Fermi level, which for the doping 
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levels of interests would be near the energy —0.8 eV of the so-called extended 
saddle-point at X. The conduction band has mostly 0-Cu anti-bonding pda- 
character (O Px - Cu dx2_y2) with the bonding partner lying 10 eV lower 
in energy. The bottom of the conduction band is seen to cross and hybridize 
with a multitude of 0-Cu pdn-hands lying below -1.2 eV. The top of the 
conduction band hybridizes strongly with a broad 0-Ca bonding pdir (O Px 
-Ca dxy) band near A. In this situation, one clearly does not want to use the 
rather ill-defined and very long-ranged Wannier orbital for describing the low- 
energy electronic structure. Rather, one wants an orbital which describes the 
band (including its dependence on other relevant low-energy excitations such 
as spin- fluctuations and phonons) in the ±200meV range around ep, that is 
an NMTO with all channels, except Cu dx2-y2 , downfolded and with as short 
a range as possible. The four dotted bands shown in each of the sub-figures 
result from such calculations ||2^. In all cases, the screening-sphere radius of 
Cu dx2-y2 was taken to be 0.62i. The upper figures illustrate a problem with 
the 3rd-generation tangent LMTO method: If Si, is taken where we want it 
to be, at the —0.8 eV saddle-point deep down in the anti-bonding pda-hand, 
then the method develops a schizophrenia near the top of the band, above 1 
eV and near M, which is apparently sufficiently far away from that the 
LMTO 'might consider' describing the bonding rather than the anti-bonding 
state. 

The resulting orbital has very long range due to the high Fourier components 
caused by the schizophrenia and, as a result, we are forced to take at a 
higher energy than we actually want. With e,j=-0.3 eV, we still get long 
range as seen in the upper left-hand figure, and in order to cure that problem 
we need to go to £,^=+0.3 eV, but then the description of the bottom of the 
anti-bonding band, the extended saddle-point in particular, has substantially 
deteriorated. In the lower left-hand figure we have now switched from the 
tangent to the chord LMTO, and that is seen to help considerably. Finally, 
the lower right-hand figure presents what might be called an 'overkill': We 
have used the discrete CMTO (iV=3) method, and the agreement with the 
exact result is superb. 



Using integrals of divided differences of MTOs. In all previous deriva- 
tions of the variational LMTO method, the LMTO was expressed as a matrix 
Taylor series and the Hamiltonian and overlap matrices (|^) were worked 

out using expressions ( pl| ) for (^(j) \ (pj and ^0 | 0^ . 

The same may be done for the general, discrete NMTO method, although 
the number of terms in the resulting series increases quadratically with N. 
For this, we first use a divided-difference form -such as (p7|)- for the NMTO 
and then need expressions for the overlap integrals, (</> [0..A^] | (p [0..M]) , and 
Hamiltonians, {(f>[O..N]\?i\ 4>[0--M]) , between divided differences of kinked 
partial waves. Since expressions (|6^) and (p^) are formally equivalent, we 
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Fig. 11. Conduction band of CaCu02 calculated by massive downfolding to a 
single Cu — NMTO (dotted) compared with the full band structure (solid) 
0. See text. 

find that, analogous to (|9l|), 

(0[O..M] I 0[O...iV]) = {(j)[Q...N]\(j)[Q..M]) = K[[Q..M].N] (101) 

,m (N) [N) (M) (M+N+l) 

\M! ' iV!/ \iV!'Af!/ (M + iV+1)!' 

where we have assumed M < N. From this result for M — N, it follows that 
the odd-ordered Hermite divided differences of the kink matrix are positive 
definite. For a contracted mesh, this overlap matrix is seen to depend only 
on M + N. 



Developing the MTO Formalism 57 



For the matrix elements of the Hamiltonian we must use: 

(0 [0..M] \H - e„| 4> [0...N]) = [0..M] \ (t> [0..n - 1, n + \..N]) 

_ {K [[0..n - 1, n + 1.. min (M, N)] n.. max (M, TV)] 
~ \ k [[0.. min (A//, TV)] ..n-l,n+ 1.. max (TW, TV)] 

(M) (TV) (M) (W-l) (M+W) 



(102) 



TW! ' TV! / \ Af! ' (TV- 1)! / (Af + TV)!' 

where the upper and lower results on the second line correspond to n ~ 
min(Af, TV). Here again, for a condensed mesh the Hamiltonian matrix de- 
pends only on Tkf + TV. 

The resulting expressions for (x^^-* | X^^'') (x^^'' 1^ ~ £n \ X^^'') con- 
tain the above-mentioned integrals times products of — em-i)- 
matrices. These expressions are by far not as explicit as equations ( p3|) and 
(M), and they are more complicated for a discrete than for a condensed mesh. 



We shall now consider a more useful application of (101)- (102) 



Charge density and total energy: Si phase diagram. The wave func- 
tion obtained from a variational calculation is: (r) = x Ci , where we 
have dropped the superscript (TV) on the NMTO. The eigen(column)vector, 
Ci, of the generalized eigenvalue equation (|^) should be normalized accord- 
ing to: cj (x I x) Ci' — Shi , or -regarding CRL^i as a matrix- according to: 
(x I x) c = 1. The charge density is now given by (p7|), which to a very 
good approximation is ( |63| ) with the energy-dependent wave functions in ex- 
pressions (|64[)-(|65|) substituted by their matrix Lagrange or Newton series. 
The computer code would use the Lagrange form: 

p (r) = X (r) cc'I'x (rY = ^ (t)n (r) L,,cc^ L\, 0„/ (r)''' , 

nn' 

SO that in this case, the density-of-states matrix F [e) in ( [66| ) should be 
substituted by: 

(occ \ 
Y^CcUlI. (103) 

Equations (|64|)-(|65|) then become: 

P^{V) = Y,Y.Y. (^«) rRL,n;R'L',n' ^R'L',n' (rji')* , (104) 

RR' LL' nn' 

Pr = X! ^L' (f) X! "^-R''" rRL,n:RL',n' ^Rl' ,n' (r) , 

LL' n7i' 

Pr (r) =T.Yl (f) Yl, (f) ^ (r) rflL,„;ijL..„. ^?j,,,„, (r) . 
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If one feels that, with the variational NMTO method, the KKR equations 
have been solved with sufficient accuracy, then one may even use (p4[)-(p6|) 
as they stand, and interpolate the energy dependences of the wave functions 
using the classical Lagrange or Newton methods ( |145D and (14(;). 

In order to solve Poisson's equation and to compute the Coulomb- and 
exchange-correlation integrals for the total energy and forces, we need to fit 
the charge density by suitable functions. The properties of p{r) to which 
we have most easy access are its spherical-harmonics expansions around the 
various sites. For the fitting we therefore choose atom-centered NMTO-like 
functions which have the following advantages: (1) they are the unitary func- 
tions for continuous fitting at non-touching a-spheres, (2) they are localized, 
(3) we know the result of operating on them with V^, and (4) the integral of 
any p r oduc t of two such functions is the energy derivative of a kink matrix 




Fig. 12. Total energy of Si as a function of the atomic volume for different struc- 
tures calculated with the full-potential LMTO method [^] and with the present 
full-charge scheme See text. 



Our fitting procedure Q can be outhned as follows: We first place a 
set of screening spheres around each atomic site. This defines our screened 
Hankel functions ( ]2^ ) and divides space into non overlapping intra-sphere 
parts and an interstitial part. It is not necessary to place screening spheres 
at interstitial sites, even though the resulting interstitial can be very large. 
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In the intra-sphere region we use a spherical-harmonics expansion of the 
charge density, with the components pRL(r) known on a radial mesh. As 
the screening spheres are relatively small this summation can be truncated 
at Z=3 or 4. In the interstitial we expand in the screened Hankel functions, 
n'^j^ (e, tr) , normalized as in (Q) and with 3 different, negative energies, of 
which the lowest is about 4 times the work function, that is: 

2 

~J2J2^RL ^r) ^RL;n = ^ TI^rl PRL + (105) 

Ti=0 RL RL 

RL (i"fl) PRL (a) + ([01] , tr) ^ Xrl,r'L' Pr'L' 

R'L' J 

for all vr > URL- With three energies, we can in principle fit continuously 
with continuous 1st and 2nd derivatives. However, in practice it is difficult 
to compute the 2nd radial derivatives of the high-Z components of the charge 
density. We therefore determine the matrix X in such a way that the fitting is 
continuous and once differentiable, that is: X = i?° [01] ^ {d{p{a)} — Bq) . 
The functions n^^ (r^) in ( |105| ) are those linear combinations of the three 
(e„, r)'s whose value and radial slope vanish in all channels at the screening 
spheres. These functions therefore peak in the interstitial region and their 
coefficients prl are determined by a least squar es fi t in the region interstitial 




t o th e MT-spheres, by sampling the full charge ( |104| ), as well as the expansion 



( 105 ) at a set of pseudo-random points. Once the expansion is obtained, 
it is very easy to solve Poisson's equation. In the intra-sphere part this is 
done numerically and in the interstitial analytically by virtue of the screened 
Hankel functions solving the wave equation. The same expansion procedure 
can be applied to the exchange-correlation energy density e(r) and potential 
/x(r). This gives a full potential. The total energy Etot is also easy to evaluate. 
The interstitial part of the integrals reduces simply to a summation over 
Hermite divided differences of the slope matrix. 

We have applied this procedure to look at the total energy of various 
possible structures for silicon ]^ . For each structure we perform a standard 
self consistent LMTO-ASA calculation. In the last iteration an expansion of 
the full charge density is made and Etot evaluated correctly. The result is 
shown in Fig. [l^ w here, for comparison, we show the full-potential LMTO 
result from RefT"^ . 



Overlapping MT-potential: Si without empty spheres. The phase dia- 
gram of Si just shown was calculated using LMTOs defined for MT-potentials 
with empty spheres. We now consider the possibility offered by Eq. ( p7| ) of 
allowing the atom-centered sphere a substantial overlap -like the 50% radial 
overlap shown in Figs. |^^- and, hence, of getting rid of the empty spheres. 

The first question is: How to construct such a potential? Our answer 
is M] that the potential should be constructed such as to minimize the 




mean squared deviation of the valence-band energies from the ones for the 
full potential. From this condition, it then follows that the overlapping MT- 
potential, '^^v (r^j) , should be the least-squares approximation to the full- 
potential, V (r) , weighted with the valence charge density. This yields a set 
of coupled equations for the shape, f (r) = v (r) — g, and the zero, g, of the 
MT-potential. The equation which arises from requiring stationarity with 
respect to 5g is of course: / {V — 'Y^v) pcPr = 0, and it means that the error 
in the sum of the valence-band energies should vanish to leading order. The 
other equations, which arise by requiring stationarity with respect to 5f (r) , 
are coupled integral equations, which are complicated due to the presence 
of the charge-density weighting. Taking the charge density to be constant in 
space, corresponds to minimizing the mean squared energy-deviation for the 
entire spectrum, rather than merely for the valence band. Now, in our present 
implementation, we only took the spatial behavior of the charge density into 
account in the 5(7-equation. The resulting potentials for diamond-structured 
Si were shown in Figs. 10 and 11 of Ref. We have recently succeeded in 
obtaining the overlapping MT potential from the full potential obtained from 
the charge density ( |105D but in the present paper wc shall only show 
results obtained by taking the full potential to be the Si+E ASA potential 
-like in Ref. 
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L r X 

Fig. 14. Band structure of Si calculated with the 3rd-generation LMTO method 
for the self-consistent Si+E MT-potential (dashed) and for the Si-centered, 60%- 
overlapping MT-approximation to it (solid). The latter calculation included the 
correction for the kinetic-energy error Eq. q27\ ) in the LMTO Hamiltonian, and the 
value of the MT-zero was adjusted in such a way that the average energy of the 
valence band was correct. Hence, the solid band structure corresponds to the last 
point on the curve marked 'ideal' in Fig. ^ |42|^ . See text. 

Fig. |l^ shows three different results for the rms error of the valence-band 
energy as a function of the linear overlap, uj = (s/t) — 1. For the overlap 
increasing up to about 30%, the rms error falls in all cases, simply because 
the overlapping MT-potential becomes an increasingly better approximation 
to the full potential. Without any overlap correction, the kinetic-energy error 
(p7|), which is of second order in the potential overlap, initially rises pro- 
portional to V (s)'^ Lu^ pO| , and this is seen to limit the maximum overlap 
to about 30%. We may, however, use the LMTO equivalent Q of Eq. 
to correct each band energy, Si (k) , and the results are shown by the two 
other curves. The dashed curve -marked 'present technique'- uses the Sg- 
equation as given above, whereas the 'ideal' curve was obtained by adjusting 
g -iteratively, because g enters the df (r) equations- to have the mean error 
of the valence-band energy vanish exactly. It is possible to improve upon the 
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'present technique' without knowing the valence-band energy a priori, and 
we are currently including charge-density weighting in the 5f (r)-equations. 
This makes the curve flatten out -like the one marked 'ideal' p3| . 

The solid curves in Fig. ^ show the Si band structure obtained with 
the 60% overlapping MT-potential, including the LMTO overlap correction, 
and determining g to yield vanishing mean error of the valence band. The 
dotted curve is the 'exact' result as obtained with a (3rd-generation) LMTO 
calculation for the Si+E potential. The errors seen in the valence band are 
certainly no larger than 30 meV, but those in the conduction band are larger. 



4 Energy-dependent linear transformations 

If one considers Fig. |l|, it might seem as if the energy-window over which an 
NMTO set yields good approximations to the wave functions will be wider if 
one starts out from energy-dependent linear combinations of kinked partial 
waves: 

0(e,r) = 0(e,r)f (e), (106) 

which have smoother energy dependencies. Normalized kinked partial waves 
and Lowdin orthonormalized kinked partial waves are examples of cases 
where the divergences of the kinked partial waves at the energies, e^^, where 
a node passes through the screening radius, are avoided. The transformation 
given by the - in general non-Hermitian ~ matrix T (e) mixes kinked par- 
tial waves with the same energy and different RUs linearly. Although the 
Hilbert spaces spanned by the energy-dependent sets, (j) {s, r) and (e, r) , are 
identical, it is not obvious that those spanned by the respective polynomial 
approximations, x*'''^'' (i") X^^^ (r) , are also identical, particularly not if 
one bears only Fig. Q in mind. 

Depending on the transformation, the resulting (j) (e, r) may completely 
have lost its original i?L-character. Since the linear combination, (j) {e, r) , of 
kinked partial waves has active radial functions on other sites, as well as at 
its own site for other L's, it is not a kinked partial wave in the usual sense, 
that is, one which could have been obtained by a screening transformation. 
Remember, that for 3rd-generation kinked partial waves, a screening trans- 
formation is not linear. In the following, we shall assume that the screening 
radii have been chosen at the previous step, in the screening calculation for 
the structure matrix, and perhaps by subsequent re-screening of the G^'s 
using (H). 

A further motivation for considering transformed kinked partial waves is 
that they might provide the freedom to obtain energy matrices (^) which 
are Hermitian. This would simplify the finite-difference expressions (8^) and 



|7|) for the NMTO so that they take the simpler form (0) which then -like 
in (^- could be diagonaHzed to leading order by the eigenvectors of E^'^\ 
From expression (p3) for the transfer matrix, we realize that the condition 
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that a transformed i?^^^^ be a Hamiltonian matrix, is that we can find a 
transformation with the property that 

>(M) I ^(Af-1)\ ^ 



^l-J = 1. (107) 

This formahsm could therefore also be the basis for obtaining an orthonormal 
NMTO set. 

Let us finally express the important equations (^6|)-(^2|) in terms of the 
transformed kinked partial waves: 

(H-e)^ (e, r) = -5 (r) K (e) , {n~e)4> (e, r) G (e) -6 (r) , (108) 

where we have defined the non-Hermitian matrices 

K{e) = K{e)f{e), G (e) = K {ey^ = f {ey^ G (e) . (109) 

Note that these definitions do not correspond to similarity transformations. 
The kink matrix, K (e) , and thereby its inverse, G (e) , were originally defined 
in such a way that they are Hermitian, but they are inherently 'skew', because 
(108) tells us that it is the 'one-sided' contraction of the Green function, 

7 (e, r) = (e, r) G (e) = 4> {e, r) G (e) , (110) 

which is invariant. For the same reason, the integrals of the products of two 
contracted Green functions, with possibly different energies, form an overlap 
matrix, 

G{s)^{He) I ke'))G{e')=- ^^'lZ^}''\ (111) 

which is independent of T (e) . 

Adding to the discussion following (^ij) about the meaning of the matrix 
equation (0„ | <pn') = | (f>n) , note that this equation does not hold in a 
general representation: <^0„ | 0„,^ = ^Tl^^I"^ (^4>n' \ (i>r}jf-^fn' ^ ^0„' | 

unless Tn = Tn' . But it is of course always true that | 0„' ^ — (4>n' \ 4>nJ ■ 
We now come to derive NMTOs from the transformed kinked partial 
waves (106). Since the arguments around expression ( |68| ) concerned the con- 
tracted Green function, which according to (110) is invariant, ( |6^ ) is im- 
changed but should be rewritten in the form: 

N 

X^^) (e, r) G (£) - (e, r) G (e) - ^ k (r) G„ (e) . (112) 

n=0 



As a consequence, (|69|) should be substituted by: 



-1 



,(^) Z\^0(r)G/' A^G \ _A^dp{r)G( A^G 



A[O..N] \A[O..N] A[O..N] \A[O..N] 



(113) 
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The last equation (113) shows that the polynomial approximation to the 
transformed energy-dependent NMTO, x*"^"* (^i r) = x^^'' (Sj r) T (e) , is 

xi^^ (r) = x'^^ (r) G [0...N] G [0...7V]"' , (114) 

which is a linear transformation. Hence, regardless of the energy-dependent 
transformation T (e) of the kinked partial waves, all NMTO sets span the 
same Hilbert space and all energy-windows are therefore identical. This 
disproves the above-mentioned naive conclusion drawn from Fig , p]. Since 



G (e) = T{e)G{e), we may express the NMTO transformation (114) as a 
Newton series (^) for f (e) : 

G [0...N] G [0...iV]"^ = (f g) [0...iV] G [0...7V]"^ (115) 

= G[M..N] G[O...Nr^ 

^fo + ..+f [0...N] (^(1) - EN-l) .. (^^"^^ - So) . 

Since the contracted Green function is invariant, so are equations (^) and 
( p2[ ) which relate the overlap and Haniiltonian integrals of such functions to 
Hermitc divided differences of G (e) . For the NMTO overlap and Hamiltonian 
matrices, we therefore obtain (|9^) and (p^), with the prefactor substituted 
by G [O..A'']~^^ , the postfactor substituted by G [O..A^]~^ , and the Hcrmite 
divided differences of G (e) unaltered. 



The first equation (113) shows that the expressions derived previously for 
the NMTOs, excluding those for integrals over NMTOs, may be taken over, 
after these expressions have been subject to the following substitutions: 

0(£,r)->0(£,r), K{e)^K{e), li''^ , . 
X (£, r) ^ X (e, r) , G [e) ^ G [e) , E^^^) ^ ^(*^) . ^ ' 

Remember, that the substitutions for K (e) and G (e) do not correspond to 
a similarity transformation. 

As long as we only consider T (e)-transformations which are independent 
of iV, the step-down relation ( [79| ) holds for the transformed NMTOs and for 
its transfer matrices, because the derivation merely made use of (|5^), which 



transforms into (108). This shows that E'^^^ — £q equals —Kq = —KqTq, as 
expected, but that: | x^"^^) = 1 does not hold. The hatted version of 

(11) therefore only holds for iV > 1. For iV = : 



\n - £o| X^°^) = -flKo% = fl (^(") - £o) ^ - £0. (117) 



The expressions for the transformed NMTO in terms of divided differences 
of transformed kinked partial waves are the hatted versions of ( |85| ) and (p7|). 
One should remember that the divided difference, 0QO..M] , r) , is a linear 
combination of the Ad+l functions 0o (r) To, .., <j)M (r) Tm, and hence, a linear 
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combination of the M + 1 divided differences: (r) , .., ([0..M] , r). This is 

the generahzation of the property: dcj) (e, r) T (e) /de\^ = (r) T + (r) T, 
used in the 2nd-generation LMTO formahsm. Exphcitly: 

([0...M] , r) = f; , (118) 

M 

= ^ {[m..M] , r) f [0..m] = (/) ([0...M] , r) Tq + .. + 0m (r) f [0...M] . 

m=0 



The transformed versions of the resuhs (101), (102) are comphcated, unless 
T (e) is independent of e. In that case, the right-hand sides just have K (e) 
substituted by f^K {e)f = K (e) . 

Usually ^0 [0..M] | [0..iV]^ ^ {^0 [0..N] \ ^ [O..M]j , unless f (e) = f , or 

the matrix is diagonal; ^(/)[O..Af] | <j>[O..N]^ = (-^[O-.A^] | (^[O-.A/J^ of course 
always holds. 



5 Hamiltonian energy matrices and orthonormal sets 



Having seen that an energy-dependent, linear transformation (106) of the 
MTO set does not change the Hilbert space spanned by the set of energy- 
independent NMTOs, but merely the individual basis functions, we now turn 
to the objective of finding a representation in which the energy matrices E^-^-^"^ 
-but not necessarily the Green matrix G (e)- are Hermitian. The energy 
matrices will then be the two-center Hamiltonians entering expressions like 
(^ for the orbitals. From (|8^), we obviously want: 

E^''^ - SM - {x^''^ \n - SmI X^''^) ^ H^''^ - SM (119) 

for 1 < M <N, and since this condition leads to the near-orthonormality 
condition ( |l07|) , it guides the way to make one of the NMTO sets -let us call 
it the Lth- orthonormal. 

In order to solve the N near-orthonormality conditions for the Hamilto- 
nian matrices, we first insert the transformed version of expression ( p^ for 
the inverse of the Afth divided difference of the Green matrix in terms of the 
transfer matrices and _ defined by (pT7|), 



-G [0...M]-' ^ f (^(0) _ (^(1) _ (^(M) _ , (120) 

into the transformed version of expression ( p^ ) for the Hamiltonian in terms 
of the 2Ai^th Hermite divided difference of the original Green matrix G (e) . 
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We then use (119) and notice that one factor ~ eu cancels out so that 

the equation may be solved for this highest transfer matrix: 



H 



(M) 



£m — 



£m-i 



-G[[0..M-1] M])To 



-If 



£m-i 



for M > I. Solving recursively for the transfer matrices, and including (117) 
at the top, we obtain the following results: 



Hi") 


-£o = 








- £1 = 




[]0]G[[0] ir'G[[]Q]f,-'^ 


^(2) 


- £2 = 




0]-'G[[0] l]G[[01]2]-iG[[0] l]G[[]0]-'fo 


fjiM) _ 


- £m = 




ft) G[[]0]<"'^ ...G[[0..M-1]M]"' 



.G[[]0] 



^ n 



(121) 



where for reasons of systematics we have used the notation (153): 
explained in the Appendix. 



The divided differences (120) of the tra nsfor med Green matrix are needed 
for specification of the transformation via (10£), the orbitals via (114), or the 
transformed kinked partial waves via (110), and are seen to be given by: 



- G[[]0]-'fo 
G[01]'' = -G[[0] l]-'G[[]0]f„" 



G[0] 
5- [01] 
G [012] 



It 



(122) 



G[[01]2]-^G[[0] 1] G[[]0]-'fo 



G [0...A/]"' = {-r G [[0..M - 1] M]-^ ...G [[ ] 0] 



Since we originally had the N+1 matrices To-.-T^r at our disposal and have 
used N to satisfy the near-orthonormality conditions, we have one, Tq, left. 
This -and thereby implicitly also the other T„'s- may now be chosen equal 
to a matrix, Tq, which makes the Lth set orthonormal. Note that whereas 
the transformation T (e) did not depend on the order of any basis set, the 
transformation T (e) does; it depends on L. 

Let us first discuss whether the transformation ( p7| ) to an orthonormal- 
ized NMTO set may at all be arrived at by an energy-dependent linear 



transformation of the kinked partial waves: According to (114), othonor- 
mality of the Lth set happens for any transformation T (e) which satisfies: 
(f-^G) [0...L] = (-G [[O...L]])^^^ , where G [[G..L]] is the {2L + l)st Hermite 
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divided difference (151) of the original Green matrix. Hence, this is a hnear 
equation between the L + 1 values of the matrix T (e)~ at the first L + 1 
mesh points, and it is therefore plausible that it may be used to fix Tq. 

The better way of writing this equation is, like for the Hamiltonian matrix. 



- -Tf ' G [[ ] 0]^"'^" ' ' ..G [[0..L]] ..G [[ ] 0]<"'^^^' f(|"''^(^^^ 

We see that the equation (x*-^-* | x^^^*) = 1, hi contrast to the equation: 
(x^*^^ |7i — Ea/I X^*^"*) = -ff*-*^-* — Em, is quadratic in all Hamiltonians, and 
therefore can only be solved by taking the square root of a matrix. 

Hence, our strategy is to choose a Tq, which makes the non-orthonormality, 

^) I x*-^') - 1 = d'-^\ (124) 



to insert (120) for G [0..L] ^ into the transformed version of expression ( p3| 
for the overlap matrix. As a result: 

(x(^) I x("') = {h^'^ - Sl) .. {h^'^ - si) - so) X (123) 

fo-i (-G [[0..L]]) f^-it (^(0) „ (^(1) _ .. (^(i) _ 



so small, that we may use an expansion like (|9^) to find Tq and the corre- 
sponding Hamiltonians iJ^*^). Of these, if^^^ equals the variational Hamilto- 
nian (Pq) with N substituted by L, and its eigenvalues are therefore correct 



to order 2i + 1. Expression (123) now tells us that: 



-^(_l)i- + i(t)^+i / (i) ^(L)\ f,(-l)^ + ^(t)"' _ 7S(-l)"'+nt)^ + '7S(-l)^ + Ht)" 



which may be solved to yield: 

(-1)^+1 



To =ToVl + 0(i 




(125) 



Here, the upper result is for L odd and the lower for L even. Since 0(-^) wiU be 
chosen s mall, and fo r L > 1 is usually of order (e^ — ei) (e^ — £q) , as we shall 
argue in (136) and (143), this transformation preserves the i?L-character of 
each NMTO. T he H amiltonian matrix (121) is seen to transform like the 
overlap matrix (123) with M substituted for L and, as a consequence. 



jjiM) 



1 + 



(M) 



-(-1)^ 



(126) 
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Similarly, from (122): 



G[O...M]"' = G[O...A/]"' \/l + 0(^ 



-(-1)^ 



(127) 



A procedure for comp uting [1 + 0] ^ , which is more robust than the 



matrix Taylor series (125), is included in our codes 



Choosing Tq. Since the near-orthonormality conditions ( |107D merely fix the 
geometrical average \ X^'^^~^^') of successive sets, the nearly orthonor- 



mal scheme (121)-(123) only makes sense if the transformation Tq of the 
kinked partial waves at Eq is chosen in such a way that the non-orthonormality 
O*-"^ is small compared with the unit matrix. The nearly-orthonormal scheme 
alone, does not make the orthonormalization integrals (x''^^'' I X^*^'*) converge 
towards the unit matrix, but make them behave like: 



This alternates with fluctuations depending on the size of (x''^'' | x''^'') • 
The first thing to do is therefore to renormalize the MTOs in such a way 



that Tg^ \\(t)RL\ J = 1> instead of (|46|). Hence, the first choice is: 

To" = (fco)"' (128) 
where fcg is the energy-mdependent diagonal matrix with elements 

l'^flx(£o)|') = K^L.RLi^o) ^ kRL,RLi£0)- (129) 

Another choice is to start with a Lowdin orthonormalized Oth-order set: 

_ i_ 

= (fcg) 'VTTo^^' (130) 

where is the non-orthonormality of the Oth-order, renormalized MTO set: 
This choice therefore corresponds to taking L = Q. 



Test case: GaAs. We have tested this orthonormalization method for GaAs 
using the minimal Ga spd As sp basis set and going all the way up to L = 3, 
that is, to a CMTO basis with the properties that = (x^^^ |H| x^^^) and 
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{x^^^ I X^^"*) = 1> so that H^^^ is a 7th-order Hamiltonian. H^'^^ and /f'^^ are 



of lower order, however, and neither of the three Hamiltonians commute. 

We diagonahzed H^^^ for L — 1,2,3 and compared with the band struc- 
tures obtained with the corresponding non-orthonormal variational method 
discussed in Sect. 3.2. Both starting choices (128) and ( |13C| ) were tried, and 
both gave fast convergence of the square-root expansions. The first choice 
which only requires evaluation of a square root at the last stage (126) but 
whose non-orthonormality O'^^^ is larger, was found to be the fastestp3. 



Aleph-representation. The renormalization (|128| ) is of the same nature as 
-but simpler than (due to lack of energy dependence)- the one performed in 



Subsection 2.3, where we went from phase-shift normalization to screening- 
sphere normalization. That diagonal transformation was given by ( ^4[ ) for the 
screened spherical waves, by (^) and (^) for the Oth-order MTOs, and by 
( ^8| ) for the KKR matrix. Since we distinguished between those two normal- 
izations by using respectively Greek and Latin superscripts for the screening, 
e.g. a and a, and since it is irrelevant, whether one arrives at a nearly or- 
thonormal representation from quantities normalized one-or- another way, it is 



logical to label quantities having the integral normalization (128) by Hebraic 
superscripts, e.g. H as corresponding to the same screening as a and a. Al- 
though not diagonal, and therefore influencin g the shape of the kinked partial 



waves, also the Lowdin orthonormalization (13C) is an energy- mdependent 



similarity transformation, and so is any of the following transformations: 

(132) 

(e) = ^ (e) (e) ee f^-^ G« (e) Tq" ^ 

with Tq arbitrary. From the latter energy-independent similarity transforma- 
tion of G (e) , the non-Hermitian matrices i^^' and E^-^\ which are given 
in terms of G (e) by respectively ([72|) and (|80|), are seen to transform like: 



if/"^) = f^-^Ll^^'Hs and = f^-^E-'^^'^fS. (133) 



This -( |l32| )-( |l33| )- has all concerned an energy-independent similarity trans- 
formation of wn-hatted quantities. 

In order to ensure that the hatted quantities are independent of which 
representation -a or H- we start out from, e.g. 

(e, r) = (e, r) = ^ {e, r) f'^ (e) = </>^ (e, r) (e) 

and 

G^ (e) = G"^ (e) =7"^ (ey^ G"" {e)f^ [eY^^ =f^ (e)"' {e)f^ (e)"'^ 
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where, from the latter, it follows that 

fN(A') _ fa(N) , pK(M) _ pa(M) 

it suffices to satisfy the relation: 

f^{e) = f^-^f^ie) , which leads to : f^^l. (134) 

In conclusion, under the substitution a — s- H, all previous equations remain 
valid, and the factors Tq may be deleted. 

The virtue of this notation is that, once we have decided upon the nor- 
malization and the screening, we can drop the superscripts; and this is what 
we shall do: From now on, and throughout the remainder of this paper, 
un-hatted quantities, i.e. the kinked partial waves, the kink and the Green 
matrices, and the Lagrange and energy m atric es, are all supposed to have 
the integral (ortho)normalization ( |12^ ) or ( |130| ), that is, they are all in the 
Aleph-representation. All equations derived previously are then unchanged, 
and To may be dropped. 

Accuracies of Hamiltonians. The accuracies of the Hamiltonians depend 
on the sizes of the corresponding non-orthonormalities. Specifically, since the 
residual error of the one-electron energy after use of the variational principle 
(|) for the set x^*"^' (r), 

is proportional to (e^ — Eq)^ ■• (ei ~ £a/)^ , neglect of the non-orthonormality, 
leads to the error: 

5ef " = (e. - sm) a|f ' + o[ie,- sof .. (e. - £m)'} , (135) 

where O^*^'' = v\o'^^^Vi and O means at the order of. The goal should thus 
be to reduce the non-orthonormality to: 

because in that case, the error from non-orthonormality will be of the same 
order as that of the residual error. This can usually only achieved for M — L. 

The order of the non-orthonormality may be found by use of the difference 
function: 

^(M)(r)_^(M-i)(^) = ,^([01],r)(ffW- 7/(^^-1))+ <^([012],r) 
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obtained from (^) and where we should take = if m < 1. As a result: 

0(M)^^^(M)|^(M)„^(M-1)^ 



(136) 



0[O1] |0[Ol])(i?(*^-i)-ei 



^(M) 



i [012] 



which is usuaUy of order - ei^ - Eq) when M > 1. 

To evaluate integrals like ((j)o\(j) [01] ) we must transform to the original 



representation using (118) and then use (|101|). In this way we get: 



00 I [01] )^{^o\^ [01]) + (00 I 0i) T [01] - K [[0] 1]+K [01] T [01] . 

(137) 



Remember, that we are using the Aleph-normahzation (132), because this 
influences the right-hand sides. For a condensed mesh, (137) reduces to: 



We shall conclude this study of the accuracy of the Hamiltonians in Eq. ( 144 ) 
below. 



6 Connecting back to the ASA formalism 

What remains to be demonstrated is that the NMTO sets, x'^-* (r) , x*^^"* (r) , 
and X*-^-* (r) , of which the two former are based on Lowdin-orthonormalized 



kinked partial waves at the first mesh point (13C), and the last corresponds to 
the L=l-set being orthonormal, are the generalizations to overlapping MT- 
potentials, arbitrary A^, and discrete meshes of the well-known LMTO-ASA 
sets given in the Overview by respectively (|l|), (0), and (^). 

Since in the present paper we have not made use of the ASA, but merely 
a MT-potential -plus redefinition of the partial waves followed by a Lowdin- 
orthonormalization- we merely need to show that the formalism developed 
above reduces to the one given in the Overview for the case A'^^l and a con- 
densed mesh. In order to bridge the gap between the new and old formalisms, 
a bit more will be done though. 
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N — 0, L — 0. For the Oth-order set we have: 
^(0) (r) = ^(0) (r) ^ (r) ^ 0„ (r) . 

All un-hatted quantities in the present section will correspond to using kinked 
partial waves, transformed to be orthonormal at this first mesh point, Eq. 



That is: All un-hatted quantities are in the Aleph-representation (132)- (134) 
with Tq given by ( 130 ). In this representation all previously derived relations 
hold, and in addition: 

To = 1 and ka = 1. (138) 

Relating back to the Overview, this means that instead of the ASA-relation 
(p^, we have ( |132D with Tq given by ( |130| ). The latter is the proper definition 
of Kq ^^^^, now that Kq = {(j)^ \ 4>q) is no longer diagonal. We now see that 
the un-hatted quantities used in the Overview were, in fact, in the Aleph 
representation. 

The overlap and Hamiltonian matrices for the Oth-order set are thus: 



h) = U"^ I X^"^ 



1 



(x'"' \n - eol X<°^) = (x'"' \n - eo\ x(°)) = H("^ - eo = ~Ko, (139) 

and with the Oth-order set being orthonormal, the Hamiltonian is variational. 
Hence, i?^"-* = is the first-order, two-center, TB Hamiltonian of the 3rd- 
generation scheme. 



TV = 1, £ = 0. For the LMTO set we have: 

X^^) (r) = </>o (r) + <t> ([01] , r) [e<^'^ - eo) ^ </> (r) + ^ (r) 



where i?^^-' -as given by ( p9[)- i s seen to become the Hermitian, first-order 
Hamiltonian H^^'> given by (139) if the mesh condenses. This proves dl). 

The Hamiltonian and overlap matrices were given in respectively (p5) and 
(p6[), and using now K — I together with (101), we see that for a condensed 
mesh 



H 



(0) 



H 



(0) 



and 
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which are exactly (|6|). Merely 

integrals like in the ASA. 

The nearly orthonormal LMTO set is: 



is not a diagonal matrix of radial 



xW(r) = <^o(r) + <^([01],r)(ffW-eo), 
and the two conditions: (x'"^^ I X^°'') = 1 = {x'"^^ I X^^'') > therefore lead to: 
'^[01] I 0o) =0= (<^o I <^[01]), and (^i | ,^o) = 1 = 



Of these matrix equations, the first means that any (pRL ([01] , i") is orthog- 
onal to any cjiR^L' (^Oji") . As a consequence, the leading term of the non- 
orthonormality ( |136| ) vanishes. The non-orthonormality of this LMTO set is 
then: 



6^'^ = (7?(l)-£o)(<^[Ol]|0[Ol])(i?«-£o) 



(140) 



which by use of (135) shows that the errors of the iJ^^^ -eigenvalues are: 



5e. 



i(i) 



0[O1] I (^[01]) (e, - £i) (e, - eo)' 



(141) 



This is one order better than the error cx (e^ — Eq) obtained by diagonal- 
ization of H'^^\ but one order worse than the error cx (e.; — ei) (e.; — eo) 
obtained variationally using the LMTO set. Hence, ij'^^ is a second-order 
Hamiltonian. From (121): 



i?(i)-£i = -GoG[[0] l]"^Go ^ -G 



G = 



1 - K- 



2! 

which for a condensed mesh is exactly (|^) 



i-K) = [i+(h(°)-£.)(0U)]" 

id mesh is exactly (|^). 
For the transformation (114) from the x to the x-set, we get by use of 



(122) 



G[01]G[01]"^ = -G[01]G[[0] l]"^Go 



-G 



G = G' 



G 



which -since from (101) 



is exactly (|^). 
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The transformation (118) of the kinked partial waves is most easily found 
by using the orthogonality of cjiQ (r) and cj) ([01] , r) together with (137). For 
a condensed mesh, the result is simple: 



(r) = 0(r) + (/)(r)T = (r) - (r) ( | 0) = 0(r)-0(r)-|, 



and well known -see Eqs. (^ and (|0). For a discrete mesh, things look more 
complicated in i^- language: From (137), 



f [01] 



-K[01]~^ K[[0] 1] 



-K [01] 



_i l-i^[01] 



£0 — £l 

where the 2nd equation has been obtained by use of (|153|): F [[0] 1] = 



_ -Fo--F[01] 



together with: Kq = 1. For ( |118D we thus obtain: 

0([Ol],r) = 0([Ol],r) + <^i(r)f [01] 

= 0([Ol],r)(l + (£i-eo)T[01]) + 0o (r) f [01] 

= 0([Ol],r)i^[Ol]-i + 0o(r)f [01] 

= {0([Ol],r) + 0o(r)r[Ol]i^[Ol]}if [01]-^ 

= {0 ([01] , r) - ^0 (r) K [[0] 1]} K [01]-' (142) 



where from (101): K [[0] 1] = (0o ] (p [01]) is the equivalent to the usual radial 
integral and the new factor K [01] in the transformation is caused by the 
presence of (pi (r) rather than ipQ (r) on the right-hand side of the top line in 

(El)- 

In order to complete the identification of the nearly-orthonormal LMTO 
representation for a discrete mesh with the ASA version (^ and (^Tj) , we need 
an explicit expression for the third parameter, which is the matrix entering 
the non-orthonormality (140). With the help of (142), and remembering that 
(po (r) and (p ([01] , r) are orthogonal, we get: 



[01] 10 [01]) = A' [01] ^^^[01] 10[O1]; 

- K [01]-' f (0 [01] 1 cb [01]) - K [[0] 1]') K [01]-' 



= K [01]-' ( K [[01]] - K [[0] 1]' ) K [Oiy 

2 



K 
3! 



K 
2[ 



where, in the third equation, we have used (101). 
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TB parametrization For tight-binding parametrizations of many bands 
over a relatively wide energy range, it is usually important to have as few 
parameters as possible. Our experience (6^,^ for the occupied and lowest 
excited bands of semiconductors and transition metals is that the off-diagonal 

elements of (0o \ = K [01] , (0o I [01]) , and (^^ [01] | 4> [01]^ may be 
neglected. This is in the spirit of the ASA. We therefore need to tabulate 
only those few diagonal elements, together with the sing le TB matrix H'^°\ 
These quantities may then be used to construct for instance the Hamiltonian 
and overlap matrices (x^"'^'' |W ^ £i| x'-"'^') and (x'"'^' | X*'"'^'') • This is like in the 
ASA, but now, we neither need this approximation nor a condensed mesh. 

AT > 1, L — Q. The nearly-orthonormal QMTO set is: 

X(^) (r) = ^0 (r) + {0 ([01] , r) + <^ ([012] , r) (i/^) - £1) } (h^^^ - Eq) 

with the non-orthonormality: 

O(') = (x^'^ I - X'^^) = (00 I [012]) - e,) - £0) + 

+ - so) {4> [10] I 4> [01]) - + .. . 

This -together with ( |135| )- shows that the eigenvalue errors of H^'^^ are: 

fc'f ^ ~Uq\4> [012]) (£, - £2) (£, - £1) (£, - £0) , 

which means, that H^"^^ is a second-order Hamiltonian like //^^^ , but different 
from it. In general, for > 1, the leading non-orthonormality is: 

6(^) « (<^o I <^ [012]) (ij(^-i) - £i) (ijW - £o) , (143) 

as seen from ( |l36| ). This means that H'^^'^ remains a 2nd-order Hamiltonian 
when A^ > 1, and that its eigenvalue errors are: 

fe^r^ « (00 I (^[012]) (£, - £jv) (£,: - £1) (e. - £0) . (144) 

This is much inferior to the variational estimate obtainable with an NMTO 
basis. Moreover, the same result would have been obtained had we started 
out from the cheaper, renormalized scheme based on ( |12^ ). Hence, with the 
present scheme only the Hamiltonians H'^^'> with M ^ L, have eigenvalues 
which are accurate approximations to the one-electron energies. 
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N — 1, L — 1. We finally use the general procedure (124)-(127) to or- 
thonornialize the nearly-orthonormal LMTO set considered above. The small 
parameter -the non-orthonormality O'^-^^^^- is thus given by (140). 

The transformation from the nearly to the completely orthonormal set is 



obtained from (127), with L = M = 1, as: 



^(1) (r) = x^i) (r) G [01] G [01]"' = x^'^ (r) 1 + O^'^ 

which is the generalization to discrete meshes and (overlapping) MT-potentials 
of the first equation (@). The resulting, orthonormal LMTO set is: 



X«(r) = 00 (r) + 0([Ol],r)(7?W-£o) 



with the third-ovder Hamiltonian obtained from (126) with L = M = 1 as: 



H 



(1) 



£l 



i + o(i) 



^(1) 



£l 



1 + 0(1 



This is the second ASA equation 



For the transformation of the kinked partial waves, we have from (125): 



00 (r) = 00 (r) 



1 + 6(1) 



and putting all of this together, we may obtain: 

0([Ol],r) « 0([Ol],r) - 00 (r) (7/(1) - eo) (0 [01] I [01] 

which is a new result. Finally, we may check that: 

{x^'^ I X^°^) = (00 I 0o) + {h'^'^ - So) (0 [01] I 0o) = 

1 + 6(1) - (7?(i) - £0) (0 [01] I [01]) (i?(i) - £0) (00 I 0o) « 1. 

7 Outlook 

Of the new developments described above, only the use of overlapping MT- 
potentials and efficient computation of total energies and forces from TB- 
LMTO-ASA charge densities were planned. Those parts turned out to be the 
hardest and still await their completion. But on the way, we did pick up a 
number of beautiful and useful instruments. Now that we have an accordion 
for playing Schrodinger, maybe Poisson can be learned as well. 
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9 Appendix: Classical Polynomial Approximations 

Lagrange and Newton interpolation. In these interpolation schemes, a 
function / (e) is approximated by that polynomial of iVth degree, Z*-^^ (e) , 
which coincides with the function at the + 1 energies, Eq, ei, .., Etv, forming 
the mesh. The error is proportional to {e ~ Eq) {e — ei) .. {s — En) ■ 

The expression for the approximating polynomial in terms of the + 1 
values of the function, / (£„) = /„, with n = 0, 1, .., N, is: 



/W(e) = ^/„ZW(£), where (e) ^ J] (^^S) 

n— m— 0,7^n 



is the Lagrange polynomial of A^th degree. It has nodes at all mesh points, 
except at the nth, where it takes the value 1. Since Lagrange interpolation 
is exact for all functions e*^ with M < N, the Lagrange polynomials satisfy 

the sum rules: e^' = J^n^o i^nf^ '"^^ (e) Jor M = 0, .., TV. 

The same approximating polynomial may be expressed as a divided dif- 
ference -or Newton- series: 

N M-1 

^/[0,..,M] n (146) 

M=Q n=0 

= / [0] + / [0, 1] (e - eo) + .. + / [0...N] (e - en-i) {e - £i) (e - eq) , 

where the square parentheses denote divided differences as defined in the 
following table: 



£0 /o = /[0] 

^Mll,/[0,l] 



In general, that is: 

rr , -, /[to,to+ l,.,n] - f[m + l,.,n,n+l] 
/ [m, m + 1, .., n, n + IJ = , (147) 

£m £n+l 

where m < n. Note that the two energies in the denominator are those which 
refer to the mesh points not common to the two divided differences in the 
nominator. Also, note their order, which defines the sign. A divided differ- 
ence, /[0...M], is thus a linear combination of /o, /i,...,/m- The divided 



differences entering ( 146 ) are those descending along the upper string in the 
table, but other forms are possible. Besides, the order of the energies need 
not be monotonic. In fact, all divided differences of degree M + 1 involving 
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M specific mesh points are identical. This means that the order of the ar- 
guments in / [0, 1, ., M — 1, M] is irrelevant, as may be seen explicitly from 
expression ( |14§| ) below. When we have a long string of arguments, we usually 
order them after increasing mesh number, for simplicity of notation. 

We may express any divided difference, / [0..A/] , entering the Newton 



form (14(:) as a linear combination of the /n's with n < M, and thereby 
establish the relation to the Lagrange form (|145| ). To do this, we apply both 
Newton and Lagrange interpolation to a function, which we take to be that 
Mth degree polynomial, /(*^) (e) , which coincides with / (e) at the first M+1 
mesh points. This is allowed, because / [0..A/] is independent of the /„'s with 
n > M. In this way, we get the identity: 

M m-1 M 



n=0 



and taking now the highest derivative, we obtain the important relation: 

M J. 

f fO-^] - E ; _ (148) 

The inverse relation, that is the expression for /„ in terms of divided differ- 



ences for a (sub)mesh containing e„, is of course just the Newton series (146) 
evaluated at the mesh point £„. 



In order to factorize {4>G) [O...A^] in expression (69) for the NMTO, we 
shall need to express the TVth-order divided difference of a product function, 
f (e) g (e) , in terms of divided differences on the same mesh of the individual 
functions. Since the product is local in energy, we start by expressing its 
divided difference in the Lagrange form ( |148| ): 

N 

{f9)[0...N] = Y^ 

n=0 llm=0,74n V^" 

For / (e) we may choose to use the divided differences in the upper, descend- 



ing string of the table. We therefore use (146) to express /„ in terms of the 
divided differences on the (0..n)-part of the mesh and thereafter reorder the 
summations: 



N N M-1 



(/g)[0...iV] = EE /[O-^^] n 



n=0 M=0 m'=0 Y[m=0,^n 

N N T-rM-l 



-E / [o-^^]E ntr'^° , ".'". g^ 

M=0 n=0 llm=0,#n i^" ^"W 
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Since Om'^o i^n -£„') = forn < M, 



N T-rA/-l / N 



n=0 nT?x=0,#n Sm) 



N y-rM-l ( _ \ 

T-lN , _ -.9" 

N 

9n 

T-rA' / _ N 

n=M llm=M,=in\Sn Sva) 



= 5[M..iV] 



according to (148). We have thus proved the binomial formula: 

N 

ifg) [0...N] =J2f [0-^'^] 9 [M..N] , (149) 



M=0 



which expresses the A^th divided difference of a product on the (O...A^)-niesh 
as a sum of products of divided differences on respectively the (0..M)- and 
(M..Af)-parts of the mesh, with M being the only point in common. Hence, 
this formula is in terms of the divided differences descending forwards along 
the upper string for /, and the divided differences descending backwards 
along the lower string for g, but this is merely one of many possibilities. For 
the special case: g (e) = e, we get the useful result: 

(e/)[0...iV] = /[0..iV-l] + eNf[O...N]. (150) 

Since the numbering of the points is irrelevant, we could of course have singled 
out any of the iV + 1 points, not merely the last. 

Newton interpolation has the conceptual advantage over Lagrange inter- 
polation that the 1st divided differences, f [n — l,n] , are the slopes of the 
chords connecting points n ~ 1 and n, and hence approximations to the 1st 
derivatives, the 2nd divided differences, / [n — 1, n, rt + 1] , are 'local' approx- 
imations to 27 times the 2nd derivatives, and so on, as expressed by (|70|). For 
the mesh condensing onto the one energy, e^, Newton interpolation becomes 
Taylor expansion, which is of course simpler. An example of this is the bi- 
nomial expression for the A^th derivative of a product: For a discrete mesh. 



there are many alternatives to (149), but for a condensed mesh, there is only 
one expression. 



Hermite interpolation. It will turn out that the NMTO Hamiltonian and 
overlap matrices are best understood and computed using the formalism of 
Hermite interpolation. Here, one seeks the polynomial of degree M + N + 1 
which fits not only the values, /„, at the -f 1 points, but also the slopes, 
/„, at a subset of M + 1 points. We shall number the points in such a way. 
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that the M + 1 points are the first. This polynomial is: 

f{M+N+l) (-^^ ^ 

^ ^ 2 ^ 1 



M 

E 

n=0 



fn+\fn-fn\ ^ 



m=0,=tn m=M+l 
N 



n=M+l 



For those interested in why this is so, here are the arguments: The product 
of Lagrange polynomials 



M 



eHe)eHe)= n 



2 N 



n 



m=M+l 



with < n < M, is of degree M + A''. At a mesh point, £ = £„/, this product 
has value 1 when < n' = n < M, value and slope when < n' ^ n < M, 
and value when M < n' < N. Since the slope is: 



M 

E 



E Sri 



,=M+1 J 



the polynomial of degree M + N + 1 : 

/ M AT 

E E 



t 771=0, m=M+l 



with < n < M, has value 1 and slope if £ = £„. If £ = £„/ £„, it has 
value and slope when < n' < M , and value and some slope when 
M <n' < N. The polynomial of degree M + N + 1 : 

(£-£„)e) (£) e) (£), 

with < n < M, vanishes at all mesh points, has slope 1 for £ = £„, slope 
for s = En' ^ En when n' and < n' < M, and some slope when M < n' < N. 
Finally, the product: 

(£) (£) 

with M < n < N, is a. polynomial of degree M + N + 1. For e = e„/ it has 
value and slope if < n' < M, value and some slope ii M < n' ^ n < N, 
and value 1 and some slope if M < n' = n < N. 



M y ^ s 2 N 

- £7n \ 1 r £ £71 



n n 



m=U ^ ^ m=M-\-l,^n 
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What we shall really need is, like in (148), (jv,/-(-jv-|-i)! times the highest 
derivative of the polynomial ^"(^^+^+1) (g). Calculated as the coefficient to 
the highest power of e, this Hermite divided difference is: 

(M+N+l) 
j(M+N+l) 




{M + N+l)\ ^ ^ n ^ 

m=0,^n m=M+l 

N 



M N 
n=M+l Yl (£„-£„)' n i£n-e„^) 

m— m—M-\-l.^n 

= lim/[0 M + N+l] = f[[O...M]..N]. (151) 

In the last line, we have indicated that the Hermite divided difference may 
be considered as the divided difference for the folded and paired mesh: 

£0 £n+i £i £n+2 ■ ■ ■ ■ £m £m+N+i ■ ■ £n 

in the limit that the energy differences, e„ = Sn+N+i — £«, between the pairs 
tend to zero. In analogy with the notation for the divided differences, we have 
denoted the (M + iV + l)st Hermite divided difference: / [[0...M] ..iV] , which 
means that the mesh points listed inside two square parentheses have both 
/„ and /„ associated with them, whereas those listed inside only one square 
parenthesis have merely /„. Like for the divided differences, the order of the 
arguments inside a square parenthesis is irrelevant, but for long strings we 
usually choose the order of increasing n. For a condensed mesh, 

(Af+AT+l) 
f 

/[[O...Ml..iVl ^ -. (152) 

As examples of Hermite divided differences we have: 
/[[0]] = /o /[[0]l]-4^ 

(153) 

/ [[01]] - ^°7,Yl°;y^^ / [[ ] o..A^] - / [0..7V] 



In the NMTO formalism the Hermite divided difference (151) comes in 
the disguise of the following double sum (^ : 
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which may, in fact, be viewed as a divided difference (14S) -albeit in two 
dimensions- but that brings httle simpHfication. So let us prove that (15f ) 
and (f54) are identical: First of all, the /„-terms of the double sum (154) 
are those for which n — n' , and they obviously equal those of the single sum 
(151). Secondly, the /„-terms in (154) are: 



fn {^n ^n' ) ^ In' (^ri 



N M 

\ ^ \ ^ In ^n' ) i In' \^n' ^ 



N r I- _ 



n=M+l Y\m=0,^n (^n ^™) n'=0 11^1=0,7^)1' (^n' 
In ^n' ) 



(155) 



M , M 

M J. N 



fn \ ^ (£"?! £^n') 



Now, according to ( |14S| ), 

M 1 

E "7"' : - -— [0...M] (156) 

is the Mth divided difference of the single-pole function 1/ (£„ — e) , provided 
that n is not on the mesh 0...M. For the sum where n is on the mesh -but 
the ri'=n-term is excluded- we have: 

M 1 M -1 

e^-e„' ^ \ ' (£„-£„/)^ 

= [0..n- l,n+ 1..M] . (157) 

(en - e) 

This result also hol ds if AI is named N, and is therefore relevant for both 
of the last terms in (155). We then need simpler expressions for the divided 
differences of the single- and double-pole functions. Guided by the results: 

1 1 1 1 d^'^ 1 M + 1 



for the derivatives, we postulate that for a discrete mesh, 

nm=o(e»-e™) (e^-e) nm=o 

(158) 



Developing the MTO Formalism 83 



For M—0, these expressions obviously reduce to the correct results, (e^ — Eq) 
and {si — £o)~ • For M > 0, our conjectures inserted on the right-hand side 



of (147) and subsequent use of ( |148| ) yield: 
-^[O..M-l]--^[l..M] 



So — 



Y\m=0 i^i ^r, 



Si ■ 



[0...M], 



[O..Af - 1] 



[1..M] 



M 



1 



which are obviously correct too. Hence, equations (15S) have been proved. 

Using finally (158) in (15(:) and (157), and right back in (155), leads to 
the /n-terms in (151). We have therefore demonstrated that: 



N AI 

EE 

n=0 n'=0 



f[n,n'] 



N 

n ( 

m—O.^n 



M 



= /[[0...M]..7V] 



(159) 



^m) W {^n' ^m') 



The final expression needed for the NMTO formalism, is one for the Her- 
mite divided difference of the product-function ef (e) . For this we can use 
(150) applied to the folded and paired mesh. As a result: 



ief)[[O...M]..N] = /[[O..M-l]..iV] + eMf[[O...M]..N]. 



(160) 



Since the numbering of the points is irrelevant, we could of course have singled 
out any of the M + 1 points, not merely the last. 
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